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CHAPTER  1 


Introduction 
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1.1  Research  Objectives  and  Methodology 

This  research  addresses  the  need  to  accurately  analyze  the  performance  of  a 
particular  class  of  frequency-modulation  (FM)  numerical  demodulation  techniques  that  are 
amenable  to  real-time  implementation.  The  class  of  demodulators  to  be  analyzed  are  those 
which  apply  the  backward  difference  estimate  of  the  time  derivative  of  the  phase  to 
produce  instantaneous  frequency  estimates.  The  receive-signal  conditions  used  for 
performance  analysis  purposes  are  that  of  a  single,  discrete-time,  constant  envelope, 
frequency-modulated  complex  exponential  in  the  presence  of  additive,  spectrally  uniform, 
Gaussian  noise.  With  the  performance  of  the  backward  difference  numerical  FM 
demodulator  as  an  established  baseline,  a  further  objective  is  to  address  the  performance 
enhancement  of  this  demodulator  in  specific  receive-signal  conditions,  by  employing 
frequency  feedback  techniques.  In  this  case,  the  receive-signal  conditions  of  interest  are 
frequency  deviations  that  are  large  relative  to  the  bandwidth  of  the  message  signal,  a  low 
sample  rate  resulting  from  decimation,  and  noise  strengths  that  yield  values  of  input 
carrier-to-noise  ratio  (CNR)  of  the  order  of  10  dB  or  more. 

The  performance  analysis  methodology  is  a  combination  of  theoretical  analyses 
and  experimental  verifications  via  non-real-time  computer  simulations.  More  correctly,  it 
is  the  noise-contaminated  complex  signal  that  is  simulated  and  used  as  input  to  the 
demodulation  processes,  for  performance  measurements  and  comparisons.  The 
demodulation  processes  themselves  are  direct  floating-point  implementations,  and  are 
more  appropriately  described  as  being  experimental  (although  non-real-time)  in  nature. 
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1.2  Background 

Since  the  pioneering  work  of  Armstrong  [1],  Carson,  and  Fry  [2],  a  wealth  of 
research  has  been  conducted  in  the  theory  and  applications  of  analog  FM  communication 
systems.  In  particular,  previous  research  has  included  various  analyses  of  the  effects  of 
noise  on  the  received  signal,  such  as  the  work  presented  by  Middleton  [3,  4],  Stumpers 
[5],  and  Rice  [6],  More  recently,  with  the  advent  of  high  speed  digital  circuitry  and 
programmable  signal  processors,  interest  has  turned  toward  discrete-time  implementations 
of  the  FM  demodulation  process. 

Various  methods  have  been  proposed  for  performing  discrete-time  FM 
demodulation.  These  include  zero-crossing  interval  averaging  for  signals  where  the 
deviation  of  the  instantaneous  frequency  from  the  carrier  is  small,  and  iterative  recovery 
techniques  using  zero-crossings  as  proposed  by  Wiley,  Schwarzlander  and  Weiner  [7] 
where  large  frequency  deviations  exist.  In  contrast  to  these  techniques,  the  research  to  be 
presented  addresses  the  problem  of  numerical  demodulation  of  uniformly  sampled 
discrete-time  signals,  using  the  complex  analytic  signal.  An  example  of  this  type  of  system 
is  the  discrete-time  FM  receiver  proposed  by  Kammeyer  [8],  in  which  sampled  in-phase 
and  quadrature  signals  are  derived  from  the  received  analog  signal  and  used  to  produce  a 
linear  estimate  of  the  original  message  modulation.  The  processing  proposed  by 
Kammeyer  is  consistent  with  the  continuous-time  definition  given  by  Cramer  and 
Leadbetter  [9]  of  instantaneous  frequency,  as  the  derivative  (with  respect  to  time)  of  the 
recovered  phase  component  of  the  analytic  signal. 
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The  earliest  work  employing  the  Hilbert  transform  for  generation  of  the  complex 
signal  is  commonly  attributed  to  Gabor  [10],  Among  the  quality  references  on  the 
generation  and  use  of  the  discrete-time  analytic  signal  are  the  works  presented  by 
Urkowitz  [11],  and  Oppenheim  and  Schafer  [12],  Consistent  with  Nuttall's  work  [13], 
Urkowitz  correctly  avoids  the  previously  held  notion  that  the  modulated  signal  must  have 
a  bandwidth  which  is  much  smaller  than  the  intermediate  frequency  (IF)  at  which 
demodulation  takes  place  within  the  receive  chain.  This  notion  is  also  quite  contrary  to 
the  common  usage  of  what  are  referred  to  as  "zero-IF"  receivers  employing  in-phase  and 
quadrature  signals,  where  such  a  bandwidth  to  IF  ratio  becomes  infinite.  Practical 
generation  of  the  sampled  analytic  signal  is  presented  by  Schussler  and  Weith  [14],  Regalia 
and  Mitra  [15],  and  Rosenkranz  [16]  using  Hilbert  transform  techniques.  Schussler  and 
Weith  present  a  design  method  for  a  causal  recursive  Hilbert  transformer.  The  work  by 
Regalia  and  Mitra  allows  for  the  use  of  known  half-band  filter  design  techniques  in 
designing  Hilbert  transform  systems.  Rosenkranz  presents  an  alternate  sampling-based 
method  to  the  more  commonplace  quadrature  downconversion  process  presented  in  [1 1], 
In  addition  to  the  reception  and  demodulation  of  man-made  FM  signals,  research  is 
also  being  conducted  on  the  concept  of  instantaneous  frequency  as  applied  to  the  analysis 
of  "naturally"  occurring  signals,  such  as  animal  vocalizations  and  human  voice.  Much  of 
this  work  also  relies  on  the  generation  of  the  sampled  analytic  signal  for  estimation  of  the 
discrete-time  instantaneous  frequency.  Among  the  more  revealing  references  for  the 
concept  of  the  instantaneous  frequency  of  a  discrete-time  signal  is  the  work  by  Boashash, 
Jones  and  O'Shea  [17],  Although  the  central  difference  is  used  rather  than  the  backward 
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difference  for  approximating  the  time  derivative,  this  work  is  an  example  where  modulo 
2 n  processing  is  applied  to  the  instantaneous  frequency  estimate  of  "monocomponent" 
signals,  taking  into  account  the  circular  nature  of  the  Fourier  transform  of  discrete-time 
signals.  (As  defined  by  the  authors,  a  monocomponent  signal  is  one  "where  there  is  only 
one  frequency  or  a  narrow  range  of  frequencies  varying  as  a  function  of  time".)  Their 
presented  applications  include  the  detection  of  harmonically  related  signals,  and 
time-varying  filtering. 

Included  in  work  since  that  performed  by  Armstrong,  Carson  and  Fry  has  been  a 
considerable  amount  of  effort  in  the  design  and  analysis  of  receive  systems  employing 
feedback  techniques.  Part  of  this  previous  research  is  the  analysis  of  devices  such  as  the 
phase-lock  loop  (PLL),  commonly  referenced  in  texts  on  communications  theory,  and  the 
less  well  known  FM  with  feedback  (FMFB)  demodulator.  The  FMFB  demodulator  as 
introduced  by  Chaffee  [18]  and  furthered  by  Enloe  [19]  and  others,  effectively  acts  as  a 
time-varying  or  tracking  filter.  However,  rather  than  changing  filter  characteristics  over 
time,  this  device  translates  the  incoming  signal  in  frequency  to  the  fixed  center  of  a 
bandpass  filter.  The  purpose  of  the  bandpass  filter  is  to  reject  out  of  band  noise,  while 
passing  the  narrow  range  of  frequencies  centered  about  the  instantaneous  frequency  of  the 
input  signal.  Inconsistencies  in  earlier  performance  comparisons  of  the  PLL  and  FMFB 
demodulators  began  to  be  reconciled  with  Develet's  research  [20],  in  which  appropriate 
conditions  lead  to  equivalent  servo-mechanisms.  Emphasis  more  recently  has  been  on  the 
discrete-time  implementation  of  these  devices. 
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The  evolution  of  the  PLL  from  continuous-time  to  various  discrete-time 
implementations  is  readily  apparent  in  the  literature.  For  example,  Rosenkranz  [21]  has 
presented  a  discrete-time  approach  to  the  PLL  for  numerical  FM  demodulation.  The 
numerical  implementation  of  a  PLL  is  often  referred  to  as  a  digital  PLL  (DPLL).  The 
demodulator  presented  by  Rosenkranz  is  a  type  of  DPLL  which  is  characterized  by  a 
non-linear  phase  detector.  Continuous-time  PLL  devices  with  a  nearly  linear  phase 
detector  characteristic  were  studied  by  Robinson  [22]  and  others,  and  were  referred  to  as 
"tanlock  loops".  As  the  name  implies,  the  characteristic  of  the  phase  detector  in  these 
devices  approached  the  piece-wise  linear  phase  characteristic  of  the  arctangent.  Research 
has  been  conducted  on  the  numerical  implementation  of  the  tanlock  loop,  as  evidenced  by 
the  work  of  Lee  and  Un  [23],  Here,  the  authors  present  a  non-uniform  sampling  device 
referred  to  as  the  digital  tanlock  loop  (DTL).  This  numerical  implementation  is  modeled 
using  an  exact  piece-wise  linear  phase  detector  characteristic.  As  a  more  recent  example, 
Ono,  Aoyama,  Hagiwara  and  Nakagawa  [24]  have  presented  a  uniform  sampling 
numerical  PLL  implementation  with  a  linear  phase  detector.  In  this  work  and  that  of  Lee 
and  Un,  we  begin  to  see  the  appropriate  trend  to  mathematically  model  the  modulo  2 n 
process,  and  include  this  process  in  subsequent  phase-domain  block  diagrams  and  models. 

Recently  there  has  been  much  interest  in  time-frequency  distribution  theory  and 
applications.  Boashash  and  White  [25]  have  employed  methods  of  instantaneous 
frequency  estimation  to  control  a  time-varying  filter.  Although  the  device  is  properly 
described  as  a  tracking  filter,  the  technique  is  more  appropriately  referred  to  as  an  open 
loop  system. 
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1.3  Overview 


As  claimed  within  the  work  presented  by  Natali  [26],  there  are  relatively  few 
publications  pertaining  to  the  automatic  frequency  control  (AFC)  loop,  a  form  of  tracking 
device  employing  frequency  feedback.  With  Natali's  research,  we  begin  to  see  an  interest 
in  the  performance  analysis  of  numerical  FM  discriminators,  as  applied  to  the  AFC  loop. 
Even  this  work,  however,  does  not  include  the  use  or  performance  analysis  of  a 
discrete-time  FM  discrimination  device  that  has  a  linear  characteristic  over  the  entire 
Nyquist  band.  Thus  it  is  evident  that  a  need  exists  to  properly  characterize  the 
performance  of  such  discrete-time  FM  demodulation  processes. 

In  view  of  the  previously  stated  objectives  of  this  research,  we  begin  in  Chapter  2 
with  the  application  of  the  concept  of  Riemann  surfaces  to  the  time-domain  phase  plane 
leading  to  the  definition  of  a  modulo  2 n  process,  referred  to  as  g[a].  Here,  a  is  a  variable 
representative  of  an  input  argument  value,  and  g[a]  is  the  output  result.  This  process  is 
necessary  for  the  proper  definition  of  the  instantaneous  frequency  of  a  modulated  complex 
exponential  sequence.  It  is  also  necessary  for  proper  theoretical  analysis  of  the 
discrete-time  angle  demodulation  processes.  This  analysis  leads  to  the  accurate  prediction 
of  changes  in  the  error  in  Riemann  sheet  identification  at  any  particular  sample  instant, 
referred  to  as  a  "phase  cycle-slip".  We  continue  in  Chapter  2  with  the  demonstration  of 
the  equivalence  of  the  traditional  quadrature  downconverter  and  the  complex  bandpass 
filter,  for  generation  of  the  in-phase  and  quadrature  sequences  from  the  original 
real-valued  input  sequence.  This  in  turn  leads  to  the  concept  of  the  "generalized 
pre-envelope",  as  presented  and  used  throughout  the  research. 
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The  concepts  developed  in  Chapter  2  are  subsequently  used  in  Chapter  3  in  the 
development  and  definition  of  the  numerical  FM  discrimination  process.  In  Chapter  3,  we 
modify  the  Arcsine  numerical  FM  demodulator  presented  by  Kammeyer  [8],  to  extend  the 
linear  phase  detector  characteristic  to  the  complete  Nyquist  frequency  band  arising  in 
complex-valued  systems.  We  also  introduce  several  other  methods  of  numerical  FM 
demodulation.  These  include  the  Arccosine,  the  Principal- Value  Arctan  (previously 
developed  by  the  author),  and  the  Standard  numerical  FM  discriminator.  The  Standard 
numerical  FM  discriminator  is  derived  from  the  conversion  to  discrete-time  of  a 
continuous-time  FM  discrimination  process.  Subsequent  demonstration  of  the  equivalence 
of  the  presented  numerical  FM  demodulators  to  the  Standard  discriminator,  and  to  the 
better  known  Direct  and  Indirect  methods,  allows  us  to  refer  to  these  devices  as  numerical 
FM  discriminators.  Chapter  3  concludes  with  the  development  of  explicit  sample  rate 
requirements  for  frequency  modulated  complex  exponential  sequences.  This  leads  to  the 
definition  of  a  pair  of  parameters  which  characterize  the  modulated  complex  exponential. 
These  parameters  are  the  relative  message  sampling  rate,  y,  and  the  discrete-time  FM 
(DTFM)  modulation  index,  8.  The  definition  of  the  process  g[a]  is  then  modified  for  use 
with  the  generalized  pre-envelope. 

In  Chapter  4,  an  analysis  is  presented  of  the  performance  of  the  backward 
difference  class  of  numerical  FM  demodulators  in  the  presence  of  additive  white  (over  the 
Nyquist  band)  Gaussian  noise  (AWGN).  Both  theoretical  and  experimental  results 
obtained  via  computer  simulations,  are  presented. 
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Finally,  in  chapter  5,  we  address  the  issue  of  performance  enhancement  of  the 
numerical  FM  discriminator  under  specific  modulation  conditions,  through  the  use  of 
specific  frequency  feedback  techniques.  Here,  the  reconstituted  numerical  FM  with 

feedback  (RNFMFB)  demodulator  is  introduced.  Enhanced  discriminator  performance 
demonstrated  via  computer  simulation  techniques. 
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CHAPTER  2 


Complex  Numbers  and  the  Generalized  Pre-envelope 


2.1  The  Geometric  Interpretation  of  Complex  Numbers 

The  purpose  of  the  proceeding  sections  is  to  review  the  mechanism  of  the  vector 
representation  of  complex  numbers  and  complex  valued  time  sequences  [27],  A 
nonrandom  complex  number,  C,  can  be  represented  mathematically  in  Cartesian  (i.e., 
rectangular)  coordinates  as 

C=Cr  +jC i  ,  (2-1) 

where  CR  is  a  real  valued  number  representative  of  the  real  constituent  of  C,  and  C}  is  a 
real  valued  number  representative  of  the  imaginary  constituent  of  C.  The  variable  j  has  the 
usual  meaning  of  the  square  root  of  -1 .  Note  that  it  is  the  association  or  pairing  of  CR  with 
Cj,  combined  with  rules  of  complex  arithmetic,  which  facilitate  the  processing  of  complex 
numbers  on  a  computing  device. 

Geometrically,  C  can  be  interpreted  as  a  point  in  a  two-dimensional  plane,  whose 
location  is  determined  by  the  ordered  pair  of  real  numbers,  {CR  ,C;},  as  depicted  in  Figure 
(2-1). 


imaginary  axis 

V 

{Q..C,} 

1  real  axis 

Or 

Figure  (2-1).  The  geometric  interpretation  of  a  nonrandom  complex  valued  constant, 
C  =  Cr+ jCi ,  as  a  point  in  two-dimensional  space. 
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It  is  also  apparent  that  this  point  in  space  is  identifiable  by  a  vector  referenced  to  the  origin 
with  a  particular  length,  L,  and  direction  angle,  a ,  as  depicted  in  Figure  (2-2) 


Figure  (2-2).  The  vector  interpretation  of  a  nonrandom  complex  valued  constant, 
C=CR+jCi ,  referenced  to  the  origin. 


The  angle  a  is  by  convention  measured  from  the  positive  real  axis.  This  angle  is 
considered  positive  when  C;  is  positive,  and  negative  when  Q  is  negative.  Also,  a  is  zero 
when  Cj  is  zero  and  CR  is  positive.  It  is  apparent  that  the  ordered  pair  of  real  numbers, 

{L  ,a  },  is  an  alternate  representation  of  the  complex  constant,  C. 

The  values  of  L  and  a  are  readily  determined  from  geometric  considerations  as 

L=JC$  +  Cj2  ,  (2-2) 


and 


a=Arctan{^  }. 


(2-3) 
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Thus  L  is  the  radial  distance  of  the  point  {CR  ,  C7 }  from  the  origin  and  is  always 
non-negative.  L  is  commonly  referred  to  as  the  envelope,  modulus  or  magnitude  of  the 
complex  constant  C,  and  is  represented  as  \C  \ .  The  angle  a  is  commonly  referred  to  as 
the  phase  or  argument  of  C ,  and  is  represented  as  ZC  .  This  phase  angle  can  be  defined 
either  as 


— 7t  <  a  <  +7t ,  i.e.,  a  e  [-7t,  +7c)  , 


or  alternatively  as 


-it  <  a  <  -hi ,  i.e.,  a  e  (-71,  -Hi]  , 

for  unambiguous  representation.  Special  note  is  made  of  Equation  (2-3),  in  which  the  first 
letter  of  the  abbreviation  for  the  arctangent  is  capitalized.  This  is  to  distinguish  this 
arctangent  of  a  ratio  of  two  known  real  quantities,  from  the  arctangent  of  a  single  known 
real  quantity.  In  the  latter  case,  an  angle  in  the  range  -n/2  to  +tc/2  results,  and  is  normally 
referred  to  as  the  principal  value.  In  the  case  of  the  arctangent  of  the  ratio  of  two 
quantities  which  are  representative  of  the  coordinates  of  a  point  in  space,  {CR ,  C, } ,  the 
signs  of  both  CR  and  C,  are  necessary  to  determine  whether  to  add  or  subtract  it  from  the 
principal  value  resulting  from  the  arctangent  of  the  ratio  of  the  coordinates.  This  in  effect 
determines  the  phase  quadrant  in  which  the  vector  lies.  Thus  the  complex  constant,  C , 
can  be  represented  in  vector  or  polar  notational  form  as  the  pairing 

C=  Iclzc. 
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As  with  the  ordered  real  pair  {CR ,  C} },  complex  arithmetic  rules  have  been  developed  for 
the  ordered  real  pair  {|C  | ,  ZC  }.  These  rules  facilitate  the  processing  of  complex 
constants  represented  in  polar  form,  on  a  computing  device. 

2.1.1  Further  Considerations  of  the  Phase  Angle  a 

From  further  geometric  considerations  it  is  evident  that  the  real  component  of  C  is 

Re{C)  =Cr=  I  Cl  cos[a]  ,  (2-4) 

and  likewise  the  imaginary  component  is 

Im{C}  =  Ci  =  ICl  sin[a]  .  (2-5) 

Employing  Euler's  identity  leads  us  to  the  exponential  representation  of  the  complex 
constant  C ,  as  the  product 


C=\C\eJa.  (2-6) 

The  exponential  form  places  in  evidence  the  fact  that  a  need  not  be  constrained  to  be  less 
than  n  in  magnitude.  This  is  a  consequence  of  the  modulo  2n  nature  of  the  cosine  and 
sine  functions  in  Equations  (2-4)  and  (2-5).  Thus  the  exponential  form  leads  to  the 
concept  of  expanding  the  complex  plane  of  Figure  (2-2)  into  a  continuum  of  "sheets" 
overlapping  in  a  helical  fashion,  which  can  be  likened  to  a  spiraling  staircase.  We  can 
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consider  each  angle  a  as  resulting  from  the  rotation  of  the  vector  C ,  starting  from  the 
phase  angle  zero.  As  the  rotation  progresses  from  slightly  less  than  +tc  to  slightly  greater 
than  +7i,  the  vector  representing  C  is  said  to  progress  from  sheet  0  to  sheet  1 .  Likewise, 
for  any  progression  across  the  negative  real  axis,  the  vector  is  considered  as  progressing 
to  the  next  sheet.  (This  helical  surface  is  a  particular  example  of  what  is  referred  to  as  a 
Riemann  Surface  [28].) 

To  allow  for  the  placement  of  the  complex  constant  on  some  particular  sheet,  r , 
the  angle  a  is  written  explicitly  in  terms  of  its  constituents  as 

a  =  ao  +  2rn  ,  r  e  integers ,  ao  e  (~7t,  +7t]  .  (2-7) 

In  so  doing,  the  various  complex  constants  represented  by  Equation  (2-6)  become  unique 
via  knowledge  of  the  integer  r .  This  integer  is  representative  of  the  number  of  angular 
rotations  of  the  vector  C  through  an  entire  2tc  cycle  as  referenced  to  the  positive  real  axis, 
and  corresponds  to  the  sheet  number  on  which  the  vector  lies.  For  example,  with  a  phase 
angle  of  -1 .  l7t ,  the  vector  C  lies  on  sheet  -1 ,  and  r  is  -1  . 

Note  that  in  many  situations,  the  integer  r  may  not  be  known.  In  this  case, 
determination  of  the  phase  angle  via  Equation  (2-3)  from  the  ordered  pair  { CR  ,Cj},  will 
yield  a0  rather  than  a  .  Thus  in  general,  the  ordered  triplet  {CR,CI,r}  is  required  for  a 
"complete"  representation  of  the  complex  constant  C .  Likewise,  the  ordered  triplet 
{|C  | ,  a0 ,  r  }  will  be  referred  to  as  a  complete  representation. 
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2.2  Processing  Considerations  in  the  Utilization  of  Complex  Numbers 


In  many  situations  where  complex  numbers  are  represented  and  processed  on  a 
computing  device,  the  representation  is  either  intentionally  or  necessarily  incomplete  in  the 
sense  that  the  integer  r  is  not  used  or  not  known.  For  example,  in  representing  the 
discrete-time  counterpart  of  a  continuous-time  bandpass  signal,  the  initial  value  of  r  is 
either  unknown  or  assigned  the  value  zero.  In  fact  if  r ,  as  a  result  of  processing,  takes  on 
another  value  other  than  zero,  it  is  often  re-assigned  the  value  zero.  In  such  a  case,  the 
process  is  utilizing  a0  rather  than  a  .  This  prompts  us  to  define  the  function,  or  more 
correctly,  the  process  g[  ]  as 

g[ a]  =  {a  +  7t  •  sgn(a)}  moduio2n  -  jc  •  sgn( a) ,  {a  +  n  •  sgn(a)} modub2*  *  0  , 

=  -Ht ,  {a  +  7t  •  sgn(a)}  moduio2K  =  0  , 

sgn(b)  =  +1  ,  b  >  0 
=  0  ,  b  =  0 

=  -l  ,b<0.  (2-8) 

A  plot  of  g[a  ]  versus  a  is  shown  in  Figure  (2-3).  The  process  g[ a  ]  is  simply 
representative  of  the  information  loss  in  going  from  the  ordered  pair  (|C  | ,  a  }  ,  triplet 
(|C  | ,  a0 ,  r  }  or  triplet  {CR ,  C7 ,  r }  representation  of  C ,  to  the  ordered  pair  (|C  | ,  a0 } 
or  {CR  ,  Cj }.  Likewise,  it  is  representative  of  the  modulo  2n  nature  of  the  quantity  e)a . 
Note  for  example  that 


aQ=g[ao\=g[cL] , 


(2-9) 
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demonstrating  that  all  phase  angles  are  mapped  to  sheet  0  ,  by  the  process  g[ a  ]. 


Figure  (2-3).  The  process  g[a  ] 

From  Equation  (2-9)  it  is  also  readily  apparent  that 

eja o  _  e;g[«o+2nc] 

_  ol  (2-10) 

Thus  for  the  incompletely  specified  complex  number  C  represented  as  {CR ,  C7 } ,  the 
equivalent  exponential  form  of  C  is  | C  '  • 

2.2.1  Implicit  Versus  Explicit  Addition  of  Phase  Angles 

We  are  now  able  to  make  a  distinction  between  the  implicit  and  the  explicit 
addition  of  the  phase  angles  of  complex  numbers.  The  addition  of  two  phase  values  which 
were  obtained  via  independent  applications  of  Equation  (2-3),  is  considered  to  be  an 
explicit  addition.  The  result  of  this  addition  ranges  from  -2tc  to  +2n  radians,  since  each  of 
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the  phases  being  summed  ranges  from  -n  to  +n  .  More  generally,  if  M  such  phase  angles 
from  sheet  0  are  added  together,  the  result  ranges  from  -M71  to  +Mjt  radians.  Note  that 
explicit  phase  angle  addition  results  in  a  phase  angle  that  is  associated  with  a  resulting 
complex  number  which  may  require  a  complete  representation;  i.e.,  the  integer  r 
associated  with  the  resultant  complex  number  may  be  non-zero. 

An  alternate  method  of  implementing  the  addition  of  the  phase  angles  of  a  set  of 
complex  numbers  is  to  do  so  implicitly,  via  complex  multiplication.  For  example,  consider 
the  resultant  product,  C ,  of  the  complex  numbers  A  and  B,  where 

A  =  \A\  g/g[ao«+2raJt] 


and 


B  =  \B\  ejg[a-ob+2rbK] 


Here,  the  phase  angles  a0a  and  a0b  and  the  integers  ra  and  rb  are  such  that 


aoa  =  g[oLa ]  and  ao*  =  g[ai]  ,  with  aa  =  oto a  +  2/*a7t  and  ab  =  aoa  +  2rbTi  . 


The  product  is  formed  as 


C-A  B 

=  {Ar  ,Aj  }  •  {Br  ,  Bi  } 

=  {Ar  •  Br  -Aj  -Bi  ,Aj- Br  +Ar  •  Bi  } 

=  {Cr,Cj}.  (2-11) 
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The  resulting  phase,  a0c ,  is  found  using  Equation  (2-3).  Equivalently  we  can  write 


C  =  A  B 

=  U|eM“‘]  •  \B\elg^ 

=  \A\- 

=  \C\e^  .  (2-12) 


Note  that  even  if  ra  and  rb  are  known,  rc  is  not  identifiable  using  the  implicit  method  of 
phase  additions.  Note  also  that  if  phase  differencing  is  desired,  the  difference  can  be 
performed  implicitly  by  first  forming  the  conjugate  of  B  as 


B*  =  { Br  ,  -Bj }  , 

and  subsequently  forming  the  product  A-B*  . 

If  it  is  desired  to  obtain  the  implicit  result  using  explicit  phase  additions,  the 
process  g[ a  ]  can  be  applied  to  the  explicit  result.  It  will  become  evident  as  this  research 
progresses  that  the  distinction  between  implicit  and  explicit  phase  additions  is  needed  to  be 
able  to  properly  analyze  the  numerical  FM  discrimination  process. 
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2.3  The  Discrete-Time  Bandpass  Signal  Representation  and  Generation  of  the 
Corresponding  Complex  Envelope 

We  turn  our  attention  now  to  the  representation  and  processing  of  signal 

waveforms  that  are  essentially  bandlimited  in  frequency,  such  that  these  continuous-time 
signals  can  be  uniformly  sampled  in  time  with  tolerable  destructive  aliasing.  For  samples 
equally  spaced  in  time  (i.e.,  uniform  sampling,)  the  Nyquist  sample  rate  requirement 
implies  that  a  lowpass  signal  of  bandwidth  IF  Hz  can  be  sampled  at  time  intervals  of 
At  =  Ts<  1/(2  W)  seconds,  to  prevent  (destructive)  aliasing  of  frequency  constituents.  This 
translates  to  selecting  a  sample  rate,  Fs  =  l/Ts ,  greater  than  twice  the  lowpass  signal 
bandwidth,  W.  For  bandpass  signals  of  approximate  bandwidth  B'  Hz,  and  originally 
centered  at  some  carrier  frequency  fc,  Hz,  sampling  can  be  performed  as  above  by  first 
downconverting  to  a  new  center  frequency  of  fc  >  B'H  Hz,  with  W=fc  +  BV2.  In  this 
manner,  we  have  not  violated  the  requirement  that  the  bandpass  bandwidth  relative  to  fc 
must  be  less  than  200%  [11],  Further  details  of  this  process  are  deferred  to  Chapter  3. 

By  sampling  such  a  bandpass,  continuous-time  signal  s(t )  at  integer  multiples  of 
At  =TS,  we  obtain  the  discrete-time  sequence  s(nT)  ,  where  the  index  n  is  an  integer 
which  identifies  a  particular  sample  of  the  newly  obtained  sequence. 

2.3.1  A  Method  of  Complex  Envelope  Generation:  The  Numerical  Quadrature 
Downconversion  Process 

The  numerical  or  digital  quadrature  downconversion  process  [1 1]  is  presented  as  a 
model  of  the  processing  which  is  required  to  generate  the  in-phase  ( i )  data 
corresponding  to  the  real  constituent  of  the  complex  envelope  sequence,  and  the 
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quadrature  (q)  data  corresponding  to  the  imaginary  constituent  of  the  complex  envelope 
sequence,  from  the  real  sequence,  s(nTs).  As  will  be  seen  in  Section  3.3,  a  properly 
downconverted  and  sampled  angle  modulated  signal  can  be  considered  to  be  a 
discrete-time  bandpass  process.  We  represent  this  input  signal ,  s(nTs ) ,  as 

s(nTs)  =  a(nTs)cos[2K/cnTs  +  <\>(nTs)]  ,  a(nTs)  >  0  ,  (2-13) 

where  a(nTs)  is  the  time  varying  amplitude  of  s(nTs) ,  and  §(nTs)  represents  the  phase  or 
angle  modulation  imposed  on  the  signal.  In  this  research,  fc  is  defined  as  the  true  center 
frequency  after  analog  downconversion  (immediately  prior  to  sampling),  which  is 
considered  to  be  close  to  the  center  of  the  Nyquist  bandwidth.  Thus  for  a  sampling  rate  of 
Fs  samples  per  second,  fc  is  about  half  way  between  0  and  FJ2  Hz,  i.e.,  fc  =  FJ4  Hz. 

As  seen  in  Figure  (2-4),  the  input  signal,  s(nT),  is  replicated  and  mixed  with  two  local 
oscillator  frequencies  which  are  in  phase  quadrature  relative  to  each  other,  and  at  the 
tuning  frequency,  f  .  This  tuning  frequency  is  set  to  exactly  (within  the  numerical 
accuracy  of  the  processing  device)  one-quarter  of  the  sampling  frequency,  since  this  is  our 
best  estimate  of  the  true  center  frequency,  fc .  The  mixing  (numerical  multiplication) 
process  generates  in  the  i  and  q  channels,  frequency  components  centered  around  both 
(/,+/)  and  (f-f)-  The  role  of  the  (identical)  lowpass  filters  is  to  remove  the 
unwanted  components  centered  around  (fc+ f).  Further  insight  into  the  workings  of  the 
numerical  quadrature  downconverter  is  arrived  at  by  analyzing  the  channels 
simultaneously  using  complex  signal  representation. 
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In-phase 


->  KnTs) 


■>  q{nTs) 


Figure  (2-4).  The  numerical  quadrature  downconverter. 

The  rules  of  complex  arithmetic  justify  the  interpretation  of  the  overall  mixing 
process  as  the  multiplication  of  the  real  input  signal  by  the  complex  exponential 


Ter*****  =  2  cos[a/0(«7’j)]  -j.  2  sin[a/0(«r,)]  ,  (2-14) 


where 


aio(nTs)  =  g[2nfi  nTs +Q]  .  (2-15) 

The  result  is  the  circular  (i.e.,  periodic)  convolution  of  the  periodic  spectrum  of  the  input 
signal  with  the  periodic  spectrum  of  the  complex  exponential  of  Equation  (2-14).  In 
effect,  the  original  signal  spectrum,  S(ei2nfT‘) ,  is  circularly  rotated  down  in  frequency  by 
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the  amount  f, .  This  is  a  direct  consequence  of  the  discrete-time  modulation  or  windowing 
theorem  [12], 

Figure  (2-5)  shows  a  representative  set  of  resulting  magnitude  spectra,  at  various 
stages  in  the  numerical  quadrature  downconversion  process.  The  pair  of  digital  or 
discrete-time  lowpass  filters  can  be  interpreted  as  a  single  real  filter  operating  on  complex 
data,  with  the  desired  magnitude  response  as  seen  in  the  figure.  Thus  after  multiplication 
by  the  complex  exponential,  but  prior  to  the  filtering,  we  have  the  pre-filtered  complex 
signal 


Sp(nTs)  —  ipinl  s)  , 


(2-16) 


where 


ipinTs)  =  2  •  s(nTs)cos[2%f i  nTs  +  0]  (2-17) 


and 


qp(nTs)  =  -2  •  s(nTs)sm[2Tifi  nTs  +  0]  .  (2-18) 

After  application  of  trigonometric  identities  and  lowpass  filtering,  we  obtain  the 
approximation  to  the  sampled  complex  envelope,  s(nTs )  ,  as 

s(nTs)  =  i(nTs)  +jq(nTs )  ,  (2-19) 


where 
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KnTs)  =  a(nTs)cos[2nfenTs  +  <$>(nTs)  -  0] 


(2-20) 


and 

q(nTs)  =  a(nTs)sm[2%fenTs  +  <k(nTs)-Q]  .  (2-21) 


Figure  (2-5).  A  set  of  representative  signal  magnitude  spectrum  plots:  a)  Representative 
input  signal  magnitude  spectrum;  b)  The  frequency  shifted  and  pre-filtered  signal 
magnitude  spectrum;  c)  The  magnitude  spectrum  of  the  resulting  complex  envelope 
estimate. 
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Here,  fe  is  the  error  in  the  estimate  of  the  true  center  frequency,  i.e.,  fe  -  (fc  -  f,)-  This 
error  is  a  consequence  of  the  noncoherent  nature  of  the  quadrature  downconversion 
process,  and  can  be  either  negative  or  positive.  In  the  absence  of  a  tuning  error,  the 
remaining  error  accounted  for  by  the  above  approximations,  is  due  to  the  non-ideal  nature 
of  the  lowpass  filters.  The  quadrature  downconverter  therefore  generates  the  sequence  of 
ordered  pairs,  {i(nTs  ),q(nTs )}  as  a  time-ordered  series  of  Cartesian  coordinates,  available 

for  subsequent  processes. 

2.3.2  A  Method  of  Complex  Envelope  Generation:  The  Complex  Bandpass  Filter 

The  processing  of  the  sequence  s(nTs )  by  the  numerical  quadrature  downconverter 
as  presented  in  Figure  (2-4)  and  described  in  Section  2.3.1,  will  be  shown  to  be  equivalent 
to  linearly  convolving  s(nT^)  with  a  complex  bandpass  filter  with  impulse  response 

heqinTs)  =  hi{nTs )  +jhq(nTs)  ,  (2-22) 

and  subsequently  (circularly)  rotating  the  result  in  frequency  by  the  amount  -ft  Hz.  To 
determine  this  impulse  response,  we  first  consider  the  complex  output,  i(nT^+jq(nT^) . 
From  Figure  (2-4),  it  is  readily  apparent  that 

i(nTs)  + jq(nTs)  =  2e~*  ■  { s(nTs )  ■  e~^> nTs }  *  hip(nTs)  .  (2-23) 

(Note  that  for  the  present  discussion,  to  simplify  analysis  we  will  use  the  standard 
approach  to  the  representation  of  the  complex  exponential,  by  not  explicitly  identifying  the 
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process  g[  ]  in  our  equations. )  Here,  hlp(nTs)  is  the  real-valued  impulse  response  of  the 

lowpass  filters  of  the  quadrature  downconverter.  Indicating  spectral  periodic  convolution 
as  ®,  the  right-hand  side  of  Equation  (2-23)  has  the  discrete-time  Fourier  transform 

2e~jd  {S(eM+W')}  •  Hip{e W')  =  2e^  {S(e^’)  •  JL  W+. f‘  +  mF*)  > 

-p-  <fi  <  y  ,  we  integers.  (2-24) 

Taking  the  inverse  Fourier  transform  of  the  right-hand  side  of  Equation  (2-24)  and 
applying  Equation  (2-23)  we  find  that 

i(nTs)  + jq(nTs)  =  2<r'®  {s(nTs)  *  [hlp(nTs)  •  eM‘nT'] }  •  e~^‘nT’  .  (2-25) 

Thus  Equations  (2-23)  and  (2-25)  demonstrate  that  the  quadrature  downconversion 
process  of  Figure  (2-4)  is  equivalent  to  convolving  s(nTs )  with  a  complex  impulse 
response,  heq(nT) ,  and  subsequently  multiplying  the  result  by  the  complex  sequence, 


2e~^e~)2l^'nTs 


The  complex  impulse  response  is  identified  as 
heq(nTs)  =  hip(nTs )  •  eflvf'nT* 

=  htp(nTs)  ■  cos[2nf inTs]  + jhip(nTs)  •  sin[27t :finTs]  .  (2-26) 

From  Equations  (2-22)  and  (2-26)  we  have 
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hj(nTs)  =  hip(nTs)cos[2nfinTs] , 


(2-27) 


and 

hq{nTs)  =  /z;^(«7’,)sin[27t//«7’,]  .  (2-28) 

This  establishes  the  impulse  response  of  the  complex  bandpass  filter.  We  now  turn  our 
attention  toward  identifying  the  complex  bandpass  filter  coefficients. 

Consider  the  lowpass  filter  impulse  response,  hlp(nT ')  ,  and  its  z-transform,  Hlp(z). 
For  Hlp(z)  of  the  form 

M 

L  bk  ■  z'k 

Hip{z)  =  -  ,  (2-29) 

1-Z  Ok  ■  Z  k 

bz  1 

the  impulse  response,  hlp(nTs) ) ,  is  determined  by  the  M+l+N  real  coefficients  bk , 
k=0,l,2,...,M  and  ak ,  k=l,2,3,...,N.  The  z-transform  of  h(nTs)  is  by  definition 

H(z)  =  S  h(nTs)  ■  z~n  .  (2-30) 

v  '  /J=-oo 

Therefore,  the  complex  sequence, 

hip{nTs)e+M'nT’ , 

has  the  z-transform 


28 


hip(nTs)e+J2*f‘ nT*  <~>  Htp(z  ■  e~^>T’ )  =  Heq(z) . 


(2-31) 


Thus  for  filters  of  the  form  presented  in  Equation  (2-29)  we  have 


tbk 

Heq(z)  =  - 

l-£  at  {z  e~j:lnflTs}-k 


2  { bk  e'MF »}** 

k=0 _ 

1-2  {ake^ikU}-z-k 


k=  1 


(2-32) 


Equation  (2-32)  implies  that  any  filter  with  the  response  as  indicated  in  Equation  (2-29) 
can  be  circularly  rotated  in  frequency  by  the  amount  f, ,  by  multiplying  the  finite  length  real 
coefficient  sequences  ak  and  bk  by  the  complex  exponential  sequence, 


Q+flnfi  kTs 


Thus  we  can  transform  the  lowpass  filter  of  the  quadrature  downconverter  into  a  complex 
bandpass  filter,  and  subsequently  linearly  convolve  the  input  sequence  with  this  complex 
filter.  The  result  of  this  convolution  can  then  be  circularly  rotated  down  in  frequency  as 
implied  in  Equation  (2-25),  to  yield  the  complex  envelope  identified  in  Section  2.3.1  . 

This  process  is  depicted  in  Figure  (2-6). 
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Figure  (2-6).  A  set  of  signal  magnitude  spectrum  plots  which  are  representative  of  the 
complex  bandpass  filter  method  of  complex  envelope  generation:  a)  Representative  input 
signal  magnitude  spectrum;  b)  The  result  of  filtering  this  input  signal  with  the  complex 
bandpass  filter;  c)  The  magnitude  spectrum  of  the  resulting  complex  envelope  estimate. 


2.4  The  Generalized  Pre-envelope 

In  light  of  the  equivalence  of  the  numerical  quadrature  downconversion  method 
described  in  Section  2.3.1  and  the  complex  bandpass  filter  /  downconversion  method 
described  in  Section  2.3.2,  the  generalized  pre-envelope  is  introduced.  Normally  the  term 
"analytic  signal"  and  the  notation  s+(nTs)  is  reserved  for  the  pre-envelope , 
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s(nTs)-e+ftltf°nT° 


to  identify  the  complex  signal  prior  to  downconversion  [29].  This  is  the  signal  whose 
spectrum  is  depicted  in  Figure  (2-6  b).  In  this  research,  however,  it  will  become  evident 
that  the  discrete-time  FM  demodulation  methods  to  be  described  need  not  rely  on  whether 
downconversion  is  explicitly  implemented  as  in  the  previously  described  methods. 
Therefore,  it  is  convenient  to  define  the  generalized  pre-envelope  as 

s+(nT„  /„)  =  ”T- ,  f  <f,  <  f- .  (2-33) 

The  pre-envelope  and  the  complex  envelope  are  simply  the  special  cases  where/,  =  / , 
and  f0-  0  respectively.  Thus  the  generalized  pre-envelope  sequence,  is  simply  the 
complex  envelope  sequence  circularly  rotated  in  frequency  to  an  arbitrary  frequency,  f0 . 
For  the  case  where  an  unknown  tuning  error,/,  exists  as  previously  described,  for 
notational  simplicity  we  will  adopt  the  convention 

s+(nTs )  =  s+(nTs,  fe )  • 

Note  that  the  real  phase  sequence  extracted  from  the  complex  generalized  pre-envelope  is 

oto  (nTsJo)  =  g[a(nTsJ0)\ ,  (2-34) 

where 
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a  (nTsJo)  =  2t if0nTs  +  §{nTs)  -  0  . 


(2-35) 


The  demodulation  methods  to  be  described  in  Chapter  3  require  as  input  the 
complex  envelope,  as  generated  by  the  preceding  techniques  or  their  equivalents.  More 
generally  it  will  be  seen  that  the  generalized  pre-envelope  can  also  be  processed  to 
estimate  the  original  message. 
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CHAPTER  3 

Analysis  of  a  Class  of  Numerical  FM  Discrimination  Methods 
Amenable  to  Real-Time  Implementation 
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3.1  The  Ideal  Discrete-Time  FM  (DTFM)  Signal 

In  Chapter  2,  methods  were  introduced  for  generating  the  discrete-time 
generalized  pre-envelope,  s+(nTs),  from  the  discrete-time  bandpass  signal,  s(nTs).  It  is 
readily  apparent  that  these  techniques  are  approximate  or  non-ideal  in  light  of  the  fact  that 
the  required  ideal  lowpass  (or  complex  bandpass)  filter  cannot  be  realized  [12].  In  fact 
there  are  numerous  distortions  of  the  original  message  and  the  modulated  signal,  which 
have  taken  place  in  all  stages  of  the  transmit/receive  process.  It  is  therefore  not  the  intent 
of  this  research  to  address  the  various  issues  regarding  the  generation  of  the  pre-envelope 
with  regard  to  its  accuracy  in  representing  the  intended  modulated  signal.  However,  it  is 
the  intent  of  this  research  to  establish  a  baseline  numerical  FM  demodulation  method,  to 
which  other  methods  can  be  compared.  To  accomplish  this  objective,  the  analyses  and 
simulations  throughout  this  research  will  utilize  what  is  considered  to  be  the  ideal 
generalized  pre-envelope, 


s+(nTs)  =  i(nTs)  +jq(nTs )  ,  (3-1) 


where 


i(nTs)  =  a(nTs)cos[2%fenTs  +  <\>(nTs)  -  0]  (3-2) 


and 


q(nTs)  =  a(nTs)sm[2nfenTs  +  $(nTs)  -  0]  .  (3-3) 
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The  class  of  numerical  FM  discrimination  methods  to  be  described  estimate  the  original 
message  by  processing  the  time  sequence  of  ordered  pairs,  {i(nTs),q(nTs)}  .  Thus 
subsequent  numerical  discrimination  methods  can  be  fairly  compared  since  all  will  be 
processing  this  same  ideal  generalized  pre-envelope. 

The  envelope,  \a(nTs)\ ,  for  the  present  will  not  be  restricted  to  being  constant,  but 
we  will  impose  the  constraint  that  a(nT)  be  greater  than  0,  as  indicated  in  Equation 
(2-13).  As  before,  fe  is  indicative  of  a  frequency  tuning  error.  The  sampled  phase  angle 
term,  <t>(w7^),  is  modeled  as  being  the  result  of  the  frequency  modulation  process 

t=nTs 

<K«r.)  =  <K0l^7-,  =  27cA/  1  m{t)dt ,  (3-4) 

o 

where  m(t)  is  the  original  bandlimited  message,  that  has  been  amplitude  limited  such  that 

-1  <  m(t)  <  +1  , 

and  is  also  zero  mean  in  time.  The  complex  signal  defined  by  Equations  (3-1)  through 
(3-4)  is  representative  of  what  will  be  referred  to  as  the  ideal  discrete-time  FM  (DTFM) 
signal.  The  constant  phase  term,  0  ,  is  comprised  of  two  sources  of  constant  phase  offset. 
The  first  is  the  arbitrary  but  unknown  phase  offset  between  the  carrier  and  the  local 
oscillator  frequencies  of  the  numerical  downconversion  process.  The  second  is  the 
constant  phase  resulting  from  the  integration  of  the  modulating  signal,  including  time  t 
less  than  0.  Both  sources  of  constant  phase  offset  are  taken  into  account  via  the  arbitrary 
but  constant  phase,  0  . 
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The  discrete-time  true  instantaneous  frequency,  y(nTs),  is  defined  as 


y{nTs)  =  fy'(nTs)  =  $'(/) I  ^nTs  =  27tA/-  m(t) I  t=„Ts  ,  (3-5) 

in  units  of  radians  per  second.  This  is  obtained  by  direct  differentiation  of  Equation  (3-4). 
(The  prime  in  Equation  (3-5)  indicates  differentiation  with  respect  to  time.)  Thus,  in  all 
approaches,  the  goal  is  to  estimate  the  sampled  instantaneous  frequency  [30],  which  is 
found  as  the  time  rate  of  change  of  phase  of  the  modulated  complex  signal,  s+(nTs). 

3.2  Numerical  FM  Discrimination  (Demodulation)  Techniques 

As  already  indicated,  the  numerical  FM  demodulation  techniques  to  be  described  in 
this  research  are  considered  to  be  of  the  same  class,  in  that  each  results  in  the  backward 
difference  approximation  to  the  time  derivative  of  the  phase,  as  a  method  of  recovering  the 
modulating  signal.  The  differentiation  of  a  continuous  function  of  time,  x(t),  can  be 
defined  as 


(3-6) 


The  backward  difference  approximation  to  the  discrete-time  signal,  x’(nTs),  resulting  from 
this  differentiation  is  [12] 

x'(nTs)  =  =  Fs .  {: x(nTs)  -  *([«  -  l]rs»  ,  (3-7) 
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and  is  readily  seen  to  be  the  (weighted)  difference  between  the  present  and  previous 
samples  of  the  sequence,  x(nTs).  The  weighting  allows  for  the  incorporation  of  the 
knowledge  of  the  value  of  Ts  into  the  approximation.  Note  that  this  approximation 
approaches  our  definition  of  the  actual  time  derivative  in  Equation  (3-6)  as  the  sampling 
interval,  Ts ,  decreases.  It  can  be  shown  (see  Appendix  A)  that  this  approximation  to  the 
time  derivative  has  the  magnitude  frequency  response 

\HB(e^) \  =  jr-  V2-2cos(®n)  ,  (3-8) 

which  for  small  values  of  |o)|  relative  to  the  sampling  rate,  closely  resembles  the  desired 
magnitude  response  of  the  time  derivative. 

These  numerical  demodulation  techniques  are  presented  in  a  progression  of 
development  beginning  with  what  is  considered  to  be  the  two  most  fundamental  methods 
of  this  class,  the  Direct  and  the  Indirect  or  Implicit  numerical  FM  discriminators. 
Following  these  methods,  the  Arcsine  method  is  arrived  at  by  employing  a  mathematical 
analysis  of  the  nearly  equivalent  continuous-time  analogy,  and  subsequently  compensating 
for  distortions  introduced  by  the  sampling  process.  The  Arccosine  method  is  presented  in 
a  similar  fashion.  As  will  be  seen,  the  Principal-Value  Arctan  approach  is  developed  via 
an  alternate  compensation  method  applied  to  the  first  stages  of  the  Arcsine  approach.  The 
last  technique  to  be  covered  in  this  chapter  will  be  referred  to  as  the  Standard  numerical 
discriminator,  since  it  is  developed  in  an  analogous  fashion  to  the  standard  or  conventional 
continuous-time  FM  discriminator. 
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3.2.1  The  Direct  and  Indirect  Methods  of  Numerical  FM  Discrimination 

Before  proceeding  further  with  the  descriptions  of  the  Direct  and  Indirect  methods 
of  numerical  FM  discrimination,  we  must  modify  our  definition  of  the  backward  difference 
approximation.  Since  we  are  applying  this  approximation  to  a  sequence  which  is 
representative  of  sequential  phase  angles  extracted  from  the  generalized  pre-envelope,  we 
must  take  into  account  the  modulo  2%  nature  of  these  angles.  Thus  we  define  the 
backward  difference  approximation  to  the  desired  discrete-time  signal,  a  '(nTs,fg),  as 

a '(nTs,f0)  =  =  ps  •  g[a0(nTs,fo)  -  a0([w  - 1  ]Ts,f0)] ,  (3-9) 

where  a 0(nTs,fo)  is  as  presented  in  Equation  (2-34).  Equation  (3-9)  is  the  backward 
difference  approximation  to  the  rate  of  change  of  phase  of  the  generalized  pre-envelope, 
s+(nTs,fo )  .  With  this  definition,  we  can  now  develop  the  Direct  and  Indirect  numerical 
discrimination  methods. 


3.2.1. 1  Direct  Numerical  FM  Discrimination 

Equations  (3-1)  through  (3-3)  represent  s+(nTs)  evaluated  at  f=fe  in  Cartesian 
form,  and  can  be  converted  to  exponential  form  using  Equations  (2-2)  and  (2-3),  resulting 
in 


s+(nTs)  =  Ji2(nTs)  +  q2(nTs)  ■  exp{/Arctan[^|]  }  .  (3-10) 
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From  simple  trigonometric  considerations,  it  is  readily  seen  that 


'i2(nTs)  +  q2(nTs)  =  \a(nTs)\ 


(3-11) 


represents  the  envelope  of  the  original  signal,  s(nTs),  as  previously  described  in  Chapter  2. 

A 

The  instantaneous  phase  message  estimate ,  <f>,  ,  at  time  sample  n  is  identified  as 

<bi(nTs)  =Arctan[|^y]  =  g[2nfenTs  +  (|)(«rs)  -  0]  =  a0 (nTs)  ,  (3-12) 

which  implies  that 

A 

— Jt  <  <| ),  <  +71  . 

In  other  words,  the  instantaneous  phase  message  estimate  will  lie  on  sheet  0  of  our 
Riemman  surface  described  in  Chapter  2.  In  a  straightforward  fashion,  Equation  (3-12) 
shows  that  it  is  possible  to  obtain  these  estimates  at  any  particular  sample  in  time,  by  using 
the  corresponding  in-phase  and  quadrature  samples.  The  Direct  method  of  numerical 
discrimination  is  simply  the  implementation  of  Equation  (3-9),  using  consecutive 
instantaneous  phase  message  estimates.  To  distinguish  between  the  explicitly  required 
modulo  2k  process  g[  ]  in  Equation  (3-9),  and  the  implicitly  applied  process  g[  ]  indicated 
in  Equation  (3-12),  the  Direct  method  of  discrimination  is  summarized  as 

MnTs)  =  btoTs)  =  Fs  gE[UnTs)-U[n-l]Ts)]  .  (3-13) 
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Here,  the  subscript  E  on  the  process  g[  ]  is  used  to  indicate  that  this  process  is  being 
explicitly  applied  to  the  backward  difference  of  instantaneous  phase  message  estimates. 
These  phase  message  estimates  have  been  obtained  via  consecutive  Arctan  calculations,  as 
in  Equation  (3-12). 

Either  the  modulo  2n  process,  gE[  ],  must  be  applied  explicitly  to  the  backward 
difference  of  consecutive  instantaneous  phase  message  samples,  or  a  "phase-difference 
quadrant  determination"  algorithm  is  required  as  follows: 

d{nTs)  =  UnTs)-h([n-\}Ts) 
if{  d(nTs )  >  +7r } 
d(tiTs )  =  d(nTs )  -  2n 
if{  d(nTs)  <  -k  } 

d(nTs )  =  d(nTs)  +  2n 

y(nTs)  =  Fs  ■  d(nTs)  (3-14) 

Note  that  the  above  quadrant  determination  algorithm  achieves  the  same  result  as 
Equation  (2-8)  only  for  the  difference  range 

-27 1  <  {<bi(nTs)  -  ^,(«rs)}  <  +271 , 

which  is  sufficient  based  on  the  range  of  phase  message  estimates  returned  by  Equation 
(2-3).  In  the  mathematically  equivalent  Indirect  approach,  it  will  become  evident  that  this 
phase  difference  quadrant  determination  is  not  explicitly  required. 
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3.2.1.2  Indirect  Numerical  FM  Discrimination 


The  Indirect  technique,  commonly  implemented  in  practice,  utilizes  the  fact  that 
multiplication  of  two  complex  numbers  in  exponential  form  is  accomplished  by  real 
multiplication  of  the  magnitudes,  and  addition  of  the  corresponding  phases.  Since  the  goal 
is  to  determine  the  difference  between  consecutive  phase  samples  as  an  instantaneous 
frequency  estimate,  it  is  possible  to  achieve  this  difference  by  multiplying  the  complex 
pre-envelope  phasor  (i.e.,  phase-vector)  at  time  sample  n,  by  the  conjugate  of  the  phasor 
at  time  sample  (n  - 1).  Thus  we  find  that 

s+(nTs)-s*+([n-\]Ts)  =  a(nTs)  a([n- l]rs)exp{/'-  [2KfeTs  +  <K«rs) 1]^)]}, 

(3-15) 

and  our  Indirect  instantaneous  frequency  estimate  is  found  as  the  "phase"  of  the  resultant 
vector  in  Equation  (3-15), 

y,(nTs)  =  F„  ■  g[2nfeTs  +  tfnT.)  -  <J>([/t  -  IF,)]  .  (3-16) 

In  units  of  radians/sec,  the  Indirect  method  of  frequency  modulation  estimation  is 
proportional  to  the  backward  difference  approximation  to  the  rate  of  change  of  phase,  plus 
the  bias  term  2nfe.  Rather  than  converting  to  exponential  form,  as  indicated  in  Section 
2.2. 1  the  multiplication  can  be  performed  in  Cartesian  form  using  implicit  phase  addition, 
since 
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s+(nT,)  -4([»-  l]r.)  =  [/(II  r.)  +/</(«  7'=)]  ■  [|([«  IR,)-w([»  IF.)]  =  A+jB, 

(3-17) 


where 


A  =  i(nTs )  •  /([«-  l]r,)  +  ^(«r,)  •  ?([»-  l]r,)  (3-18) 

is  the  resulting  real  term,  and 

B  =  q(nTs)  ■  /([«  -  1]^)  -  i(nTs)  ■  q([n  -  1  ]TS)  (3-19) 

is  the  resulting  imaginary  term.  Thus  we  obtain  our  frequency  modulation  estimate  as 

yi(nT s)  =  /vArctan[j]  .  (3-20) 

Equations  (3-18)  through  (3-20)  summarize  the  Indirect  numerical  FM  discrimination 
technique,  given  the  in-phase  and  quadrature  samples  sequences,  i(nTs)  and  q(nTs )  . 


3.2.2  Arcsine  Numerical  FM  Discrimination 

The  Arcsine  method  of  numerical  FM  demodulation  [8]  is  arrived  at  by  first 
considering  the  continuous  time  analogy  to  Equation  (3-12)  which  expresses  the 
instantaneous  phase  message  estimate  as 
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$/(<)  =  ardanljjl  =  g[2 nf.t  +  <K0  -  9]  . 


(3-21) 


The  estimate  of  the  rate  of  change  of  modulated  phase,  <j>  ■(/),  is  found  by  utilizing 

4[arctan{w(0}]  = - ,  (3-22) 

dt1  1  V"J  1  +u2(t)  dt  ’  v  ’ 


for  the  angle  range 


-f  <  arctan{«(0}  <  y  . 


From  Equations  (3-21)  and  (3-22)  it  is  apparent  that 


1  + 


fit) 


V  a  L 


'  (0 


'<0  -i- 


q«) 


di 

dt 


>  (!) 


i2(t)  +  q2(t) 


(3-23) 


for  the  instantaneous  phase  angle  range 

T  <arctan{||}  <f  . 

In  the  methods  of  numerical  FM  discrimination  being  described,  we  have  available  to  us 
sampled  versions  of  the  in-phase  and  quadrature  signals,  but  do  not  have  sampled  versions 
of  their  corresponding  time  derivatives.  The  Arcsine  method  applies  Equation  (3-7),  the 
backward  difference  approximation  of  this  time  derivative,  to  i(nTs)  and  q(nT)  to  yield 
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i(nTs)  •  \q{nTs)  -  q(\n  -  1]7,)}  -  q(nTs)  •  {i{nTs)  -  i([n  -  1]F)} 

i2(nTs)  +  q2(nTs ) 


gfor,)  •  /([/i  -  TO  -  am  •  g([”  ziM 

i2(nTs)  +  q2(nTs ) 


(3-24) 


which  can  be  scaled  by  F,  for  our  frequency  modulation  estimate. 

By  employing  Equations  (3-2)  and  (3-3)  in  Equation  (3-24)  we  find  that 


t/,  a(nTs)a(\n -  l]Fj){ sin CcosD- cos C sin D} 

M"Tl)  = - a^n TV) 


(3-25) 


where 


C  =  2; ifetiTs  +  <t )(nTs)  -  0  (3-26) 

and 

D  =  2ti fe[n -  1]TS  +  <|)([/i  -  l]r4)  -  0  .  (3-27) 


Application  of  a  trigonometric  identity  to  Equation  (3-25)  along  with  Equations  (3-26) 
and  (3-27)  yields 


T,  a{nTM\»-W,)MC-D) 

MnTt)  = - o^tT) 
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_  a{nTs)a({n  -  l]?;)sin[27t/er,  +  -  (})([»  -  1]?’5>] 

a2{nTs)  (3-28) 

Under  the  condition  that  a(nTs )  =  a([n  -  l]rs),  we  arrive  at 

to«r,)  =  sin[27t/en+(f)(«rs)-(l)([n-l]7’,)]  ,  (3-29) 


which  implies  the  need  for  an  arcsine  correction,  to  properly  recover  an  estimate  of  the 
modulating  waveform.  Thus  from  Equations  (3-24)  and  (3-29)  we  find  that 


yas{nTs )  =  Fs  ■  arcsin[^(«r5)] 


=  FS 


arC.ji1rg(»^)’»([»-1]^H("^)'g([»-1]U)i 

*■  i2(nTs)+q2(nTs)  ■* 


(3-30) 


can  be  used  to  estimate  the  time  rate  of  change  of  phase  for  the  phase  difference  range 


f  <  [27 lfeTs  +  4>(nTs)  -  <K[«  -  I]!7,)  ]  <  f 


In  comparing  the  Arcsine  method  as  summarized  by  Equation  (3-30),  to  the 
previously  described  Direct  and  Indirect  numerical  FM  discrimination  methods,  two  very 
basic  observations  can  be  made.  The  first  observation  is  with  regard  to  the  range  of  phase 
differences  over  which  the  message  can  be  recovered.  It  would  appear  that  we  are  limited 
to  differences  in  the  range  -rc/2  to  +n/2  exclusively,  to  be  able  to  utilize  Equation  (3-30). 
In  actuality,  we  can  use  the  sign  of  the  numerator  of  the  argument  in  Equation  (3-30)  and 
the  sign  of  Equation  (3-19)  to  decide  whether  to  add  -n,  0,  or  +n  to  the  negated  principal 
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value  returned  by  the  arcsine  (prior  to  scaling  by  the  sample  rate,  Fs ).  This  new 
phase-difference  quadrant  determination  algorithm  required  by  the  Arcsine  method  can  be 
summarized  by 


d{nTs)  =  arcsin[/2(MJ?;(wr;- 

if(  A  <  0  ) 
if(  B  <  0 ) 

d{nTs)  =  -n  -  d(nTs ) 
else 

d(nTs)  =  +n-d(nTs) 
yas(nTs )  =  Fs  •  dinTs)  , 


(3-31) 


where  A  and  B  are  defined  as  in  Equations  (3-18)  and  (3-19).  The  second  observation  is 
that  the  Arcsine  method  is  arrived  at  via  the  constraint  that  consecutive  envelope  samples 
are  essentially  equal  (i.e.,  the  signal  has  a  constant  modulus).  Note  that  in  the  event  that 
the  constant  modulus  condition  is  not  initially  present,  we  can  correct  the  envelope  by 
normalizing  the  in-phase  and  quadrature  signal  components,  thus  imposing  this  condition. 
This  would  be  done  at  the  expense  of  the  processing  requirements  of  the  extra  square-root 
required  for  normalization  as  implied  by  Equation  (3-11),  and  the  extra  divisions  which 
would  be  required. 
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3.2.3  Arccosine  Numerical  FM  Discrimination 

Given  the  existence  of  the  Arcsine  numerical  discrimination  method  and  the 
Indirect  method  which  requires  an  arctangent  calculation,  it  is  naturally  expected  that  a 
discrimination  method  exists  which  requires  an  arccosine  correction.  Consider  the  ratio 


A 

i\nT s)+q2(nT s) 


i(nTs}i([n-\]T,)+q(nTs)q([n-\]Ts) 

i2(nTs)+q2(nTs) 


(3-32) 


From  Equations  (3-2)  and  (3-3),  the  denominator  of  Equation  (3-32)  is 


i2(nTs)  +  q2{nTs)  =  a2(nTs )  . 


Similarly,  the  numerator  simplifies  to 

i(nTs)  ■  /([«-  \]Ts)  +  q(nTs)  ■  q([n-X\Ts)  =  a(n Ts)a([n -  l]7’i){cosCcosZ)-l-sinCsinZ)} 

=  a(nTs)a{[n  -  1  ]  Ts)cos(C  -  D ) 

=  a{nTs)a([n  -  1]^)  •  cos[27i feTs  +4)(«7’5)  -  <}>([/?  -  l]7’s)] 

Here,  C  and  D  are  as  defined  in  Equations  (3-26)  and  (3-27).  With  the  above 
substitutions,  the  instantaneous  frequency  estimate  is  equivalent  to 

Vcv11')  a\nTs) 
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For  the  constant  modulus  case,  we  define 


pac(nTs )  =  Fs  •  arccos[(j)'(/?rs)] 


=  Fs  ■  arccos 


A 

i\nTs)+q2{nTs) 


(3-34) 


The  arccosine  function  will  unambiguously  yield  a  value  between  0  and  n  radians.  The 
phase-difference  quadrant  determination  algorithm  that  is  required  is  therefore  simply 


d(nTs)  =  arccos 

if{£<0.} 

d(nTs)  =  - d(nTs ) 
yacipTs)  =  Fs  ■  d(nTs) , 


/2(«rs>f92(nr,). 


(3-35) 


where  A  and  B  are  defined  as  in  Equations  (3-18)  and  (3-19).  Thus,  the  Arccosine 
method  of  numerical  discrimination  is  as  summarized  by  Equation  (3-34)  and  Algorithm 
(3-35). 


3.2.4  Principal-Value  Arctan  Numerical  FM  Discrimination 

The  Principal-Value  Arctan  numerical  demodulation  method  can  be  developed  by 

partially  utilizing  the  development  of  the  Arcsine  method,  previously  described.  By 

making  the  observation  that  the  backward  difference  approximation  to  the  time  derivative 
as  defined  by  Equation  (3-7)  is  most  valid  at  time  t  =  [n  -  \\TS ,  or  half  way  between 
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consecutive  time  samples  (see  Appendix  A)  it  becomes  apparent  that  Equation  (3-24)  can 
be  modified  as 


VtinT.y 


q(nTs)  ■  i([n -  1]7^)~  i{nTs)  - q{[n-\]Ts) 


i(nT syHi[n-\]T s 


*] 


+ 


q(nTs)+q([n-l]T> 


(3-36) 


Essentially,  an  attempt  has  been  made  to  time-align  the  sampled  in-phase  and  quadrature 
signals  with  their  corresponding  time  derivative  estimates,  by  substituting 


i{/(nr,)  +  i([n-l]r,)} 


for  i(nTs),  and 


[{q(nT,)+q(ln-nT,)} 


for  q(nTs )  respectively.  As  can  be  seen,  the  numerator  in  Equation  (3-24)  remains 

unchanged,  and  the  denominator  is  an  estimate  of  the  square  of  the  envelope  at  time 
t=  [n-  \]TS  . 

By  applying  Equations  (3-2)  and  (3-3)  to  Equation  (3-36),  we  find  that 


&(»r,) 


4  sin[27 ifJs  +  UnTs)  -  <K[w  -  1]^)] 


a2(«7’J)+a2([»-l]r,) 

a(nT,)a([n-l]T,) 


+  2  cos[27i/e  Ts  +  (| ){nTs)  -  <|>([«  -  1  ]  r,)] 


(3-37) 


Again,  under  the  condition  that  a(nTs )  =  a([n  -  l]rs),  it  is  found  that 


§'t(nTs)  =  2tan[ 


2 


], 


(3-38) 
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after  application  of  a  trigonometric  identity.  The  modulating  signal  can  be  recovered  as 


ypv(nTs)  =  2 Fs  ■  arctan[j<j>f(«7;)]  ,  (3-39) 

A  i  ^ 

where  §t(nTs)  is  as  defined  in  Equation  (3-32). 

Note  that  the  arctangent  required  in  Equation  (3-39)  is  different  from  that  required 

in  Equation  (3-20)  since  the  range  of  the  obtained  sampled  phase  differences  has  been 

reduced  by  a  factor  of  1/2,  i.e., 


?  *  Wt(pTs)  <  f  , 

as  implied  by  Equation  (3-38).  Thus,  the  Principal-Value  Arctan  numerical  FM 
demodulation  method  as  summarized  by  Equations  (3-36)  and  (3-39)  requires  neither 
phase  quadrant  determination  as  in  the  Direct  and  Indirect  methods,  nor  phase-difference 
quadrant  determination  as  in  the  Direct,  Arcsine  and  Arccosine  methods.  Like  the  Arcsine 
and  Arccosine  methods,  this  numerical  demodulation  technique  is  best  applied  to  constant 
modulus  input  signals. 


3.2.5  Standard  Numerical  FM  Discrimination 

The  conventional  or  standard  analog  FM  discriminator  [30]  consists  of  a 
"slope-circuit"  followed  by  an  envelope  detector,  as  shown  in  Figure  (3-1  a).  The 
frequency  response  of  the  slope-circuit  is  such  that  a  constant-envelope  frequency 
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modulated  signal  at  the  input  results  in  an  output  signal  with  an  envelope  that  is 
proportional  to  the  original  message.  The  slope  circuit  is  designed  to  operate  at  a 
particular  intermediate  frequency,  fc, ,  as  seen  in  the  representative  response  in  Figure 
(3-1  b).  This  response  is  approximately  linear  over  the  frequency  interval 

(A*  -  A/)  ^  !/ 1  <  (/c/  +  A/)  . 

Typically,  the  constant  envelope  condition  is  imposed  on  the  received  FM  signal  via 
amplitude  limiting  circuitry,  to  remove  any  noise-/distortion-induced  envelope  variations. 
Thus  the  analog  FM  discriminator  converts  the  instantaneous  frequency  variations  into 
proportional  amplitude  variations,  for  subsequent  envelope  demodulation. 


constant-envelope, 


frequency  modulated 
signal 

slope-circuit 

. . . "s 

envelope 

modulation 

estimate 

_  _  N* 

(centered  at  /c' 

_ J _ ’  _  A  /.  \ 

detector 

* 

with  deviation  A f  ) 


Figure  (3-1).  a)  The  conventional  analog  discriminator,  b)  Representative  magnitude 
frequency  response  of  the  slope-circuit. 
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The  Standard  numerical  FM  discrimination  method  is  the  discrete-time 


implementation  of  the  conventional  analog  device.  In  the  discrete-time  case,  we  will  be 
processing  the  ideal  generalized  pre-envelope,  s+(nTs)  ,  centered  at  the  frequency  f=fe, 
i.e.,  the  slightly  mis-tuned  complex  envelope.  Of  the  various  stages  of  processing 
required,  it  is  the  discrete-time  counterpart  to  the  slope-circuit  which  requires  some 
detailed  development. 

We  borrow  upon  previous  developments  in  both  Chapter  2  and  Appendix  A,  to 
arrive  at  a  numerical  slope-circuit  implementation.  From  Figure  (A.  3 -2),  it  is  apparent 
that  the  response  of  the  backward  difference  identified  as  HB{&  Q)  approximates  that  of  an 
ideal  slope-circuit,  HA(t e* Q) ,  over  the  frequency  interval  0  to  FJ2  Hz.  It  is  also  apparent 
that  if  the  generalized  pre-envelope  were  centered  near  Fs  /4  Hz,  it  could  be  filtered  using 
the  (real)  backward  difference  impulse  response,  to  generate  a  signal  with  an  envelope 
modulation  which  is  proportional  to  the  instantaneous  frequency.  Note  that  the  allowable 
range  on  the  instantaneous  frequency  is  from  0  to  Fs  12  . 

Rather  than  rotating  the  generalized  pre-envelope  from  fe  Hz  to  (fe  +  Fs  /4)  Hz,  the 
impulse  response  of  the  backward  difference  can  be  (circularly)  rotated  in  frequency  by  the 
amount  -Fs  /4  Hz,  using  the  methods  described  in  Section  2.3.2.  The  result  is  a  (complex) 
impulse  response  which  can  be  used  to  directly  filter  the  complex  envelope,  s+(nTs).  The 
impulse  response  of  the  backward  difference  filter  is 


Ts  ■  hB(nTs )  =  +1,  n  =  0, 

=  0,  else.  (3-40) 
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Frequency  rotation  of  these  coefficients  is  accomplished  by  multiplying  the  sequence  of 
Equation  (3-40)  by  the  complex  exponential, 


e  J[2%'  *'nTs  ]  =  cos[2tc  •  ^  •  nTs ]  -j  sin[27t  •  •  nTs ] 

=  cos[|  •  n]  -7  sin[|  • «].  (3-41) 


The  result  is  the  discrete-time  slope-circuit  complex  impulse  response, 


Ts  ■  hsi(nTs)  =  +1,  n  =  0, 

=  +7,  «  =  -l, 

=  0  ,else.  (3-42) 


This  discrete-time  slope-circuit  is  implemented  in  a  straightforward  fashion  as  shown  in 
Figure  (3-2). 


Figure  (3-2).  The  discrete-time  slope-circuit  derived  from  the  backward  difference  filter. 
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The  output  of  this  slope-circuit  is 


s+(nTs)  + j  ■  s+([n  - 1  ]TS)  =  [ i(nTs )  -  q([n  -  1]^)]  +j  •  [q(pTs)  +  i([n  -  IF,)]  (3-43) 

If  we  consider  forming  the  square  of  the  envelope  of  the  signal  of  Equation  (3-43),  the 
result  is 

l5+(«r,)+y-5+([n-  l]r,)l2  =  [i\nTs)  +  q2(nTs)]  +  [i2([n-l]Ts)  +  q2([n-l]Ts)] 

-2  •  [i(nTs)  •  q([n  -  1]^)  -  q(nTs)  ■  /([»  -  \]TS)]  .  (3-44) 

As  already  indicated,  in  the  conventional  analog  discriminator,  a  constant  envelope 
condition  is  imposed  via  amplitude  limiting  circuitry.  In  a  similar  fashion,  a  normalized 
envelope  condition  can  be  imposed  on  the  complex  envelope,  prior  to  the  slope-circuit. 
Under  this  condition,  the  first  two  terms  in  brackets  in  Equation  (3-44)  are  each  equal  to 
unity.  In  this  case  Equation  (3-44)  simplifies  to 

k(«rs)  +j  •  s+([»  -  l]Ts)\2„ormalized  =  2  ■  [1  +  q(nTs)  ■  i([n  -  1  ]TS)  -  i(nTs )  •  q([n  -  1]^)]  . 

(3-45) 

By  simply  subtracting  the  constant  term  and  scaling  by  a  factor  of  1/2,  the  Standard 
numerical  FM  discrimination  estimate  becomes 

i',(nTt)  =  q(nTs)  •  i([n  -  1]T,)  -  i(nTs)  ■  q{[n  -  1]^)  (3-46) 
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Figure  (3-3)  summarizes  the  Standard  numerical  FM  discrimination  processing  steps. 

In  comparing  Equations  (3-46)  with  the  intermediate  step  in  the  Arcsine 
discrimination  method  of  Equation  (3-24),  it  is  noted  that  under  the  normalized  envelope 
condition,  these  equations  are  identical.  It  is  then  immediately  evident  that  once  again,  an 
arcsine  correction  is  required  as  indicated  in  Equation  (3-30)  to  achieve  a  linear  estimate 
of  the  original  message.  Likewise,  to  increase  the  range  of  allowable  instantaneous 
frequency  values  to  the  entire  (complex)  Nyquist  band,  the  phase-difference  quadrant 
determination  algorithm  in  (3-31)  can  be  employed. 


complex  data 


»  real  data 


Figure  (3-3).  The  Standard  numerical  FM  discrimination  method. 
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Thus,  the  Standard  numerical  FM  discrimination  method  is  found  to  be 
mathematically  equivalent  to  an  intermediate  step  in  the  Arcsine  method.  With  the 
appropriate  corrections  previously  described,  the  output  of  the  Standard  discrimination 
process  will  be  equivalent  to  the  output  of  the  Arcsine  method.  This  is  interesting  to  note, 
given  that  the  Arcsine  and  the  Standard  methods  were  developed  using  very  distinct 
approaches.  In  fact,  it  will  be  demonstrated  in  the  following  section  that  the  Direct, 
Indirect,  Arcsine,  Arccosine,  Principal- Value  Arctan  and  the  (corrected)  Standard 
methods  of  numerical  FM  demodulation  are  all  mathematically  equivalent. 


3.2.6  Discrimination  Method  Comparison 

In  Section  3.2.5  it  has  been  established  that  the  distortion-corrected  Standard 
numerical  discriminator  and  the  Arcsine  numerical  discriminator  are  mathematically 
equivalent.  We  will  subsequently  establish  the  additional  result  that  the  Arcsine  method 
(and  therefore  the  distortion-corrected  Standard  method)  is  also  equivalent  to  the  Direct, 
Indirect,  Arccosine  and  Principal-Value  Arctan  discrimination  techniques.  Thus  the 
numerical  demodulation  techniques  presented  in  Section  3.2  are  considered  to  be  of  the 
class  of  FM  demodulators  which  employ  the  backward  difference  approximation  to  the 
time  rate  of  change  of  phase,  as  the  message  estimate.  Further,  since  these  numerical 
demodulators  include  the  distortion-corrected  Standard  discriminator,  justification  of  the 
term  "discriminator"  for  each  member  of  this  class  is  established.  Each  member  is 
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therefore  considered  to  be  of  the  class  of  backward  difference  numerical  FM 
discriminators. 

We  first  consider  the  equivalence  of  the  Direct  and  Indirect  processing  techniques. 
Equivalence  can  be  established  by  demonstrating  that  the  instantaneous  frequencies 
specified  by  Equations  (3-13)  and  (3-16)  are  the  same.  From  Equation  (3-12),  this  implies 
that  we  must  demonstrate  the  equivalence  of  the  Direct  process, 

yd{nTs)  =  Fs  •  gE[g[2ifenTs  +  $(nTs) - 0]  - g[2nfe[n  -l]Ts  +  4>([» -  IF)  -  0]]  ,  (3-47) 

and  the  Indirect  process, 

yi(nTs)  =  Fsg[2%feTs+<b(nTs)-<b([n-l]Ts)]  (3-48) 

Applying  modulo  arithmetic  rules  it  can  be  shown  that 

g[ai  -a2]  =gfe[ai] -g[a2]] ,  (3-49) 

which  immediately  demonstrates  the  equivalence  of  Equations  (3-47)  and  (3-48),  with  a, 
identified  as 


ai  =  2TtfenTs  +  §(nTs)  -  0  , 


and  a2  identified  as 


a2  =  2ti fe[n  - 1  ]TS  +  <)>([«  -  IF,)  -  0  . 
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Thus  we  find  that 


yd(nTs)  =  Fs  ■  g[2nfeTs  +  <Ki nTs )  -  <K«7^>]  .  (3-50) 

Consider  now  the  Arcsine  method  as  summarized  by  Equation  (3-30)  and  the 
phase-difference  quadrant  determination  of  Algorithm  (3-31).  From  Equations  (3-29)  and 
(3-30)  along  with  Algorithm  (3-31)  we  have 

yas(nTs)  =  Fs  g[2nfeTs  +  ^(nTs)-^[rt-l]Ts)]  .  (3-51) 

The  process  g[  ]  indicated  in  Equation  (3-51)  is  a  result  of  the  phase-difference  algorithm. 
This  algorithm  is  simply  a  variant  of  the  phase  quadrant  determination  process  required  by 
the  Arctan  introduced  in  Equation  (2-3)  of  Section  2.1.  In  comparing  Equations  (3-51) 
and  (3-48),  it  is  apparent  that  the  Arcsine  method  of  discrimination  is  mathematically 
equivalent  to  the  Indirect  (and  therefore  also  to  the  Direct)  method  of  instantaneous 
frequency  estimation.  Likewise  for  the  Arccosine  method,  Equations  (3-33),  (3-34)  and 
Algorithm  (3-35)  demonstrate  that 

yac(nTs)  =  Fs-g[2TifeTs  +  ^nTs)-^[n-l]Ts)]  .  (3-52) 

Finally,  from  the  properties  of  the  arctangent  required  in  Equation  (3-35)  and  the  results 
of  Equation  (3-34),  the  Principal-Value  Arctan  method  yields 
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ypv(nTs)  =  Fs  ■  g[2nfeTs  +  tfnT.)  - $>([n  -  1  ]TS)]  . 


(3-53) 


Thus  the  Arccosine  and  Principal-Value  Arctan  methods  are  also  equivalent  to  each  of  the 
aforementioned  techniques. 

It  has  been  demonstrated  that  the  distortion-corrected  Standard  method,  the  Direct 
and  Indirect  methods,  the  Arcsine  and  Arccosine  methods  and  the  Principal- Value  Arctan 
methods  are  all  mathematically  equivalent  numerical  FM  demodulation  techniques.  Note 
however,  that  differences  do  exist  in  the  processing  requirements  of  each  method. 
Therefore,  one  of  the  methods  will  be  selected  based  on  implementation  considerations,  as 
an  "established  baseline"  for  subsequent  performance  analysis  comparisons. 


3.2.6.1  Selection  of  a  Baseline  Numerical  Discrimination  Technique 

Each  of  the  previously  described  numerical  FM  demodulation  techniques  require, 
as  a  minimum,  a  single  division  and  a  trigonometric  inversion.  However,  in  view  of  the 
additional  processing  required  by  the  Standard,  Arcsine,  Arccosine  and  Principal- Value 
Arctan  discrimination  methods  to  impose  the  constant  modulus  condition,  these  methods 
will  not  be  considered  further.  It  is  acknowledged  that  depending  on  specific  application 
requirements,  any  one  of  these  methods  may  be  considered  "better"  than  the  others.  For 
example,  if  an  application  happens  to  require  an  envelope  estimate  in  addition  to  the 
instantaneous  frequency,  the  envelope  calculation  is  no  longer  considered  to  be  an  extra 
processing  burden.  In  this  case,  the  Standard,  Arcsine,  Arccosine  and  Principal- Value 
Arctan  methods  become  more  attractive  from  an  implementation  viewpoint. 
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Since,  for  the  purposes  of  this  research,  it  is  the  instantaneous  frequency  which  is 
of  primary  interest,  the  Direct  and  Indirect  methods  remain  as  baseline  candidates  due  to 
their  insensitivity  to  amplitude  variations.  Further  still,  due  to  its  conceptual  simplicity, 
the  Direct  method  is  arbitrarily  selected  as  the  chosen  baseline  numerical  FM 
discriminator.  It  is  pointed  out,  however,  that  since  each  of  these  methods  are 
mathematically  equivalent,  subsequent  performance  analyses  and  simulations  of  the  Direct 
method  apply  equally  well  to  the  remaining  discrimination  techniques. 
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3.3  Sampling  Rate  Considerations 

The  sampling  theorem  as  put  into  general  practice  in  the  field  of  communications 
engineering  by  Nyquist  [31]  and  formalized  by  Shannon  [32],  has  become  commonplace  in 
the  scientific  and  technical  community.  Simply  stated  the  sampling  theorem  indicates  that 
a  real  lowpass  signal  having  spectral  content  limited  to  the  frequency  range  OtofV  Hz, 
can  be  accurately  represented  by  instantaneous  samples  of  the  amplitude  of  the  signal 
taken  at  time  intervals  of  A t=Ts<  \/(2W)  seconds.  Entire  texts  (e.g.,  [12])  have  been 
devoted  to  the  processing  of  the  resulting  discrete-time  signal,  and  to  the  reconstruction 
into  continuous  time  of  the  extracted  information  or  control  bearing  signals.  (In  fact, 
generalized  versions  of  the  sampling  theorem  do  not  require  uniformly  time-spaced 
samples,  and  also  includes,  for  example,  sampling  the  signal  and  its  derivatives  to  form  a 
discrete-time  representation  [32][33].) 

Shannon's  presentation  of  the  sampling  theorem  also  allows  for  the  conversion  of  a 
real  bandpass  (i.e.,  bandlimited)  signal,  centered  at  f=fd ,  from  continuous  to  discrete  time 
by  employing  the  proper  frequency  downconversion  techniques  as  stated  in  Section  2.3. 
Reformulating  this  result,  a  real  bandpass  signal  of  approximate  bandwidth  B'  Hz,  and 
originally  centered  at  some  carrier  frequency /c,  Hz,  can  be  sampled  by  first 
downconverting  to  a  new  center  frequency  of  fc  >  B'l 2  Hz,  with  W=fc  +  572.  An 
idealized  version  of  this  process  is  depicted  in  Figure  (3-4). 

As  Slepian  appropriately  points  out  [34],  strictly  speaking,  truly  bandlimited 
lowpass/bandpass  signals  do  not  exist  in  reality,  since  this  would  imply  that  such  a  signal  is 
present  for  all  time.  In  light  of  this,  he  established  the  concept  of  essentially  or  really 
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bandlimited  signals,  to  distinguish  between  exact  mathematical  models  and  models  that  are 
less  restrictive  for  practical  implementation.  Urkowitz  applied  this  result  and  correctly 
concludes  that  this  approximate  bandwidth  of  a  real  bandpass  signal  must  be  less  than 
200%,  relative  to  its  center  frequency,/, .  This  center  frequency  is,  as  implied,  the  center 
or  midpoint  of  the  signal  frequency  band. 


a) 


r\ 


B’<B 


/  \ 

-fe¬ 

rn 

- / - 1 

A 

- \ - 

/ 

A,  . 

— i — 1 - 

fc 

\ 

fc '  ' 

. . ~v 

B/2 ) 

V 

C fc ■ 

-B/2) 

.  anti-alias  lowpass  filter 

•> 

-v  b 

r\  . 

- - , - 

-{2fc,  -  B/2) 

-w  -f 

JC 

|  1  V 

/  W 

Jc 

(2fc,  -  B/2) 

\S(f)\ 

d) 

h) 

- - - 1  i 

-W  f 

J  C 

■  > 

f  w 

Jc 

f 


f 


f 


f 


Figure  (3-4).  Bandpass  signal  pre-acquisition  processing,  a)  A  representative  magnitude 
spectrum  of  a  bandpass  signal  of  bandwidth  B'  <B,  centered  at/,  Hz.  b)  The  magnitude 
frequency  response  of  the  analog  local  oscillator  used  for  frequency  downconversion. 
c)  The  result  of  convolving  the  original  spectrum  with  the  local  oscillator  spectrum; 

/  =  b/2,  W=fe  +572.  d)  The  representative  spectrum  of  the  desired  signal,  s(t),  prior  to 
conversion  from  analog  to  discrete  time. 
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Thus  Shannon's  previously  presented  method  of  preparing  a  bandpass  signal  for 
conversion  from  analog  to  discrete  time,  and  Urkowitz's  requirement  as  stated  above,  are 
consistent. 

We  now  turn  our  attention  to  the  representation  of  a  bandpass  signal,  s(t),  as  the 
product  of  the  real  lowpass  signal,  a(t),  and  the  real  bandlimited  signal,  cos[2tc/cH-(|)(O-0]. 
As  a  consequence  of  Bedrosian's  product  theorem  [35],  in  order  for  the  signal,  s(t), 
modeled  as 


5(0  =  a(t)  •  cos[27t/c/  +  <K0  -  0]  (3-54) 

to  be  bandpass,  it  is  required  that 


\A(f)  1=0,  \f  l>/i, 

\A(f) I  >  0,  !/l</i,  (3-55) 

and  simultaneously  that 

|?{  cos[27t/ct  +  <K0  -  0] }  |  =  0,  \f\  <f2 , 

>0,/2<!/l</3, 

=  0,  l/l>/3,  (3-56) 

with  /2  >fx .  Here,  \A(f)\  =  |^{a(0)|  is  the  magnitude  spectrum  of  the  real  lowpass  signal, 
a(t).  (Equality  in  Equations  (3-55)  and  (3-56)  can  be  interpreted  as  in  the  practical  sense 
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developed  by  Slepian.)  Figure  (3-5)  is  representative  of  this  application  of  Bedrosian's 
product  theorem,  for  the  case  where  f2  approaches  fx ,  and  a(i)  is  real.  Note  that  in  view 
of  Slepian's  results,  this  product  theorem  can  be  viewed  both  as  a  bandpass  signal 
construction  technique  and  as  justification  of  bandpass  signal  representation,  as  applied  to 
the  signal  s{t)  resulting  from  the  pre-sampling  process  depicted  in  Figure  (3-4). 
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Figure  (3-5).  Representative  spectra  depicting  a  special  case  of  Bedrosian's  product 
theorem,  a)  The  magnitude  spectrum  of  the  real  lowpass  signal,  a(t).  b)  The  magnitude 
spectrum  of  the  approximately  bandlimited  signal,  cos[a(r)],  where  a (t)  =  2%ft+§(t)-Q  . 
c)  The  magnitude  spectrum  of  the  product,  a(*)cos[<x(f)]. 
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We  conclude  that  the  pre-sampled  signal,  s(t),  resulting  from  the  proper 
downconversion  and  anti-aliasing  filtering  of  a  bandpass  signal  centered  at /c,  Hz,  retains 
the  bandpass  representation  indicated  in  Equation  (3-54).  The  resultant  discrete-time 
signal,  s(nTs),  is  therefore  (within  an  arbitrary  constant  phase  offset) 

s(nTs)  =  s(t) I t=nT,  =  a(nTs)  ■  cos[2nfcnTs  +  (j)(«7's)]  .  (3-57) 

It  is  this  real  discrete-time  "bandpass"  signal  which  is  processed  as  described  in  Sections 
2.3.1  and  2.3.2  to  generate  the  complex-valued  generalized  pre-envelope,  s+(nTs).  Thus 
for  the  bandpass  signal,  s(t),  of  bandwidth  B  Hz,  the  required  acquisition  sample  interval 
is 


i  s  ( acquisition )  ^  2 B  SCCOnds.  (3-58) 

Note  that  in  additionally  modeling  this  signal  as  being  representative  of  the  frequency 
modulation  process  as  in  Equation  (3-4),  it  is  further  required  that  a(nTs)  >  0,  as  indicated 
in  Equation  (2-13).  This  condition  prevents  the  sign  of  the  amplitude,  a(nTs),  from 
distorting  the  instantaneous  phase.  This  amplitude  has  initially  been  modeled  as  being  time 
dependent  to  allow  us  to  demonstrate  the  insensitivity  of  the  Direct  and  Indirect  numerical 
discrimination  methods  to  these  amplitude  variations.  However,  in  adopting  this  more 
general  model,  it  has  been  necessary  to  address  the  previous  bandpass  representation 
considerations.  In  subsequent  analyses  we  will  adopt  the  more  restrictive  model  in  which 
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the  amplitude  is  a  constant,  a(nTs)  =  a  In  this  constant  amplitude  case,  we  can  closely 
approximate  the  bandwidth,  B,  using  Carson's  rule  [30], 

B  =  2(Af+fm)Uz.  (3-59) 

Here,/m  is  the  highest  significant  frequency  component  in  the  message  signal,  m(t). 

3.3.1  Sampling  Rate  Requirements  of  the  Ideal  Complex  DTFM  Signal 

Once  generated,  it  is  common  practice  to  increase  the  sampling  interval  of  the 
complex  envelope  by  a  factor  of  2.  Thus  the  sampling  interval  requirement  for  the 
complex  envelope  becomes 


T s  (decimated)  —  2  ‘  Ts  (acquisiton) 


which  from  Equation  (3-58)  becomes 


(3-60) 


T,  (decimated)  ^  g  Seconds. 


(3-61) 


A  set  of  signal  spectra  representative  of  this  decimation  or  downsampling  process  is 
shown  in  Figure  (3-6).  The  ability  to  do  so  is  simply  a  consequence  of  the  lack  of  signal 
energy  in  the  frequency  bands 
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where 


*  s  ( acquisition ) 


Therefore,  increasing  the  sample  interval  by  a  factor  of  2  by  simply  discarding  every  other 
sample,  does  not  lead  to  destructive  aliasing  of  frequency  components  in  the  resulting 
decimated  complex  envelope  [12],  (As  a  consequence  of  decimation  by  a  factor  of  2,  the 
resulting  spectrum  is  1/2  the  amplitude  prior  to  decimation,  to  be  consistent  with  the 
discrete-time  Fourier  transform  definition.) 


2  2 


Figure  (3-6).  A  set  of  signal  magnitude  spectrum  depicting  the  decimation  of  the  complex 
envelope  by  a  factor  of  2.  a)  |7{  s  (nTs) }  |  with  Fs  equal  to  the  acquisition  sample  rate, 
b)  2|7{  5(«r?)}|  with  F  equal  to  the  decimated  sample  rate. 
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Likewise  the  generalized  pre-envelope,  once  generated,  can  be  properly 
represented  by  samples  at  time  intervals  of  Ts^decimated) .  This  is  evident  since  circular 
frequency  rotation  can  be  performed  either  before  or  after  the  downsampling  process, 
resulting  in  exactly  the  same  time  sequence  of  ordered  pairs,  {i(nTs),q(nT)}  Note 
however  that  this  and  any  further  decimation  in  time  must  be  done  after  generating  the 
generalized  pre-envelope  sequence.  Thus  we  do  initially  require  a  sample  interval  as 
defined  in  Equation  (3-58),  for  the  method  of  acquisition  which  has  been  defined. 

The  issue  that  needs  to  be  addressed  is  the  determination  of  the  largest  time 
interval,  resulting  from  decimation,  that  can  be  allowed.  We  do  so  by  considering  the 
requirement  to  extract  an  estimate  of  the  original  message,  m(f).  The  message  signal  is 
modeled  as  being  bandlimited  such  that 

lM(/'ra)l  =  |^{/w(0}|  >0,  Icol  <®m  , 

=  0,  Icol  >  oom  ,  CDm  =  27 zfm  .  (3-62) 

For  each  ordered  pair  {i(nT),q(nT)},  backward  difference  numerical  FM  discrimination 
will  result  in  the  estimate 


yd(nTs)  =  i'd(nTs)  =  Fs  •  g[Ts  •  hB(nTs)  *  {2nfenTs  +  $(nTs)  -  0}]  ,  (3-63) 

as  implied  by  Equation  (3-51).  Here,  hB(nT)  is  as  defined  in  Equation  (3-40)  with  the 
response  indicated  in  Appendix  A.  (The  subscript  d  can  be  thought  of  as  representing 
either  the  Direct  method  specifically,  or  more  generally,  backward  difference  numerical 
discrimination .)  From  the  results  of  Appendix  A  and  Equation  (3-5)  we  conclude  that 
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yd(nTs)  =  Fs  •  g[2nfeTs  +  Ts  •  /^7;)  *  <^'(nTs)] 

=  ^  •  gilnfeTs  +  Ts  •  hD(nTs )  *  {27iA/-  w^r*)}]  .  (3-64) 

Here,  hD(nTs)  is  the  impulse  response  of  the  distortion  transfer  function,  HD(e? n),  due  to 
the  non-ideal  nature  of  the  backward  difference  approximation  to  the  true  time  derivative. 
Equation  (3-64)  is  a  key  result,  in  that  it  indicates  the  need  to  specify  two  distinct  sample 
time-interval  requirements  for  complex  DTFM  signals. 

The  first  requirement  is  a  result  of  applying  Nyquist's  criterion  to  the  message 
signal,  m(t).  Since  the  linear  convolution  internal  to  Equation  (3-64)  does  not  result  in 
frequency  components  outside  the  range  indicated  in  Equation  (3-62),  we  have  the 
requirement 


77  >  1 

5  (< decimated)  ^  2  f 


(3-65) 


Defining  the  ratio 


Y  = 


(3-66) 


we  have  the  equivalent  requirement 


y  >  1  . 


(3-67) 


We  will  refer  to  y  as  the  relative  message  sampling  rate.  The  requirement  identified  in 
Equation  (3-67)  ensures  that  the  discrete-time  estimation  of  the  original  message,  as 
identified  in  Equation  (3-64),  is  sampled  at  a  sufficient  rate. 

The  second  sample  rate  requirement  is  a  consequence  of  the  modulo  2%  process, 
g[  ],  in  Equation  (3-64).  We  immediately  conclude  that  this  results  in  a  mapping  of  all 
phase  differences,  to  sheet  0  of  the  phase-difference  plane.  In  order  to  unambiguously 
recover  an  estimate  of  the  original  message,  it  is  additionally  required  that 

— 7C  <  2nfeTs  +  Ts  ■  hD(nTs )  *  {2k  Af-  m(nTs)}  <  +k  .  (3-68) 

Since  both/c  and  m(nT)  can  independently  be  either  positive  or  negative,  this  translates 
into  the  requirement 


Ts  (decimated)  <  2{\fe\+¥\hD(nTs)*m(nTs)\}  » 


(3-69) 


for  all  values  of  the  index  n.  From  Appendix  A  it  is  apparent  that  the  gain  of  the  distortion 
response,  |//D(ei°)|,  does  not  exceed  unity  at  any  particular  frequency.  In  addition,  by 
design  m(nTs )  does  not  exceed  unity  in  magnitude.  Therefore,  in  most  modulation 
scenarios,  Equation  (3-69)  simplifies  to  the  requirement 

Ts  0 decimated)  <  2.{|/C|+A/'}  ‘  (3-70) 
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In  this  form  it  is  evident  that  this  second  sampling  interval  requirement  ensures  that  the 
maximum  deviation  of  the  instantaneous  frequency  from  0  Hz ,  does  not  exceed  FJ2  in 
magnitude.  Thus,  this  requirement  prevents  instantaneous  frequency  aliasing.  Defining 
the  ratio 

8  =  8e  +  8m  ,  (3-71  a) 


where 


8e  = 


2  l/el 
Fs 


(3-71  b) 


and 


§OT 


2  A  / 
Fs  ’ 


(3-71  c) 


we  equivalently  require 


8  <  1  .  (3-72) 

We  will  refer  to  8  as  the  DTFM  modulation  index,  5m  as  the  DTFM message  index,  and  8e 
as  the  DTFM  tuning-error  index.  Note  that  "overmodulation"  (i.e.,  8  >1)  translates  into 
instantaneous  frequency  aliasing. 
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The  DTFM  modulation  index ,  5,  and  the  relative  message  sampling  rate,  y,  can  be 
related  to  the  modulation  index,  P,  of  the  continuous-time  FM  signal  by  considering  the 
product 

y-S  =  t  +  T-  0-73) 

9  Jm  Jm 

The  modulation  index  of  the  continuous-time  FM  signal  is  typically  defined  as  [30] 


(3-74) 


Thus  from  Equations  (3-73),  (3-74),  and  previous  results  we  obtain 

P  =  Y-S-f  =  Y-8„  ■  (3-75) 

J  m 

The  sampling  interval  requirements  summarized  by  Equations  (3-65)  and  (3-70) 
for  proper  recovery  of  the  message  signal  from  the  complex  DTFM  signal  can  be 
combined  as 


T  ,,  .  .  „ 


.  •_  r  1 


Note  that  since  we  have  modeled  the  DTFM  signal  as  having  a  constant  amplitude,  the 
decimation  process  will  not  result  in  frequency  aliasing  in  the  envelope  sequence. 
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3.3.2  Instantaneous  Frequency  Aliasing  and  the  Generalized  Pre-envelope 

In  representing  our  instantaneous  frequency  estimate  as  in  Equation  (3-64),  it  is 
implied  that  the  input  to  the  discrimination  process  is 

s+(nTs,fe )  =  a-{  cos[2n :fenTs  +  §(nTs)  -  0]  +j  sin[2n fenTs  +  §(nTs)  -  0] } .  (3-77) 

Consider  the  case  where  we  have  as  input  the  generalized  pre-envelope 

s+(nTs,f0 )  =  a  ■  {cos[2nf0nTs  +  §(nTs)  -  0]  +j  sin[2nf 0nTs  +  §(nTs)  -  0]}, 

(3-78) 

Analogous  to  Equation  (3-64),  this  results  in  the  instantaneous  frequency  estimate 

yd(nTs,f0 )  =  Fs  ■  g[2%f0Ts  +  Ts  •  hD(nTs)  *  {2nA /•  m(nTs)}\  .  (3-79) 

Note  that  for  values  of  f0  near  +/-  FJ2,  even  the  slightest  message  signal  amplitude  can 
result  in  instantaneous  frequency  aliasing.  Rather  than  discarding  the  discrimination  result 
indicated  in  Equation  (3-79),  it  is  possible  to  correct  the  frequency  rotation  implied  by  the 
term,  2%f0Ts.  Define  the  modulo  2k  process,  g[a  ,/J,  as 


where  g\  ]  is  as  in  Equation  (2-8).  By  employing  the  property  of  modulo  arithmetic 
indicated  in  Equation  (3-49)  along  with  Equations  (3-79)  and  (3-80),  it  is  readily  apparent 
that 


yd(nTs,0)=Fs  g[Ts-yd(nTs,fo),fo]  .  (3-81) 

Thus  in  the  case  that  the  center  frequency,  fa ,  is  known  exactly,  the  result  of  the 
backward  difference  discrimination  process  can  always  be  corrected.  In  effect,  the  process 
g[  a  ,/,]  changes  the  reference  of  the  process  g[a]  from  0  radians,  to  2nf0/Fs  radians,  as 
depicted  in  Figure  (3-7). 

As  indicated  in  Section  2.3.1,  in  practice  the  true  center  frequency,/, ,  is  not 
known  exactly.  This  results  in  the  tuning  error, /e ,  which  is  normally  small  relative  to 
FJ2.  In  this  case,  the  instantaneous  frequency  estimate  becomes 

yd(nTs)  =yd(nTs,fe ) 


=  Fs  •  g[Ts  -ydinTMJo -fe] .  (3-82) 

Thus  the  sampling  rate  requirements  remain  as  in  Section  3.3.1.  The  process  indicated  in 
Equation  (3-82)  will  be  referred  to  as  a  phase-difference  "downwrapping",  since  in  effect 
both  frequency  downconversion  (or  rotation)  and  phase-difference  unwrapping  or 
quadrant  determination  are  performed  simultaneously. 

In  conclusion,  the  generalized  pre-envelope  sequence,  s+(nTs,f0),  can  be  processed 
directly  by  the  numerical  discrimination  techniques  previously  described.  The  resulting 
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sequence  is  subsequently  processed  as  in  Equation  (3-82),  to  yield  the  desired 
instantaneous  frequency  estimate.  Note  that  the  discrimination  and  the  downwrapping 
process  can  be  performed  on  a  sample  by  sample  basis.  Therefore,  the  class  of  backward 
difference  numerical  FM  discriminators  are  particularly  useful  in  real-time  applications. 
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Figure  (3-7).  The  geometric  representation  of  the  process  g[ a  ,/J  as  a  rotation  of  the 
angle  reference  axes  used  in  the  modulo  27t  process  g[a].  In  this  particular  representation, 
a  positive  value  of  f0  is  depicted. 
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CHAPTER  4 


Backward  Difference  Numerical  Discrimination 
in  the  Presence  of  Additive  White  Gaussian  Noise 
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4.1  Phase  Modulation  Recovery  and  Phase  Domain  Performance  Measures 

In  Chapter  3,  various  approaches  to  the  generation  of  the  backward  difference 
approximation  of  the  rate  of  change  of  phase  were  introduced  as  instantaneous  frequency 
estimates.  From  among  these  approaches  the  Direct  method  has  been  selected  as  a 
baseline  numerical  FM  discrimination  process  for  analysis  purposes,  although  each  of 
these  approaches  results  in  the  same  estimate.  Therefore  any  reference  made  to  the 
baseline  numerical  discrimination  process  will  refer  to  the  Direct  method  specifically,  or 
more  generally,  to  any  one  of  these  equivalent  processes. 

We  summarize  the  baseline  discriminator  using  Equations  (3-12)  and  (3-13)  as 

fanTs)  =  Fs  •  gE[UnTs)  -  ♦,([»  -  1]7^)]  (4-1) 


in  units  of  radians  per  second,  where 


A 

* 


,{nTs)  =Arctan[ 


Im{ complex  envelope) 
Re{ complex  envelope) 


J. 


(4-2) 


It  is  convenient  to  define  the  phase  message  estimate  as  the  accumulation  of  the 
instantaneous  frequency  sequence  of  Equation  (4-1).  Thus  in  units  of  radians  we  define 


kd(nTs)  =Ts-{i  <\>d(kTs)  +  <M-7;)}  (4-3) 
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as  the  baseline  discriminator  phase  message  estimate.  Here,  is  an  arbitrary  initial 

phase  value  assigned  to  the  phase  at  time  t  =  -Ts .  This  is  also  the  initial  value  assigned  to 
the  instantaneous  phase  at  time  t  =  -Ts  during  the  numerical  discrimination  process,  i.e., 


Usually  the  initial  phase  4>d(-7^)  is  assigned  the  value  zero,  such  that  it  need  not  be  retained 
for  phase  message  recovery.  Note  that  the  accumulation  process  of  Equation  (4-3)  is  an 
inversion  of  the  backward  difference  calculation  [12]  and  has  the  transfer  function 
[#B(ein)]'1.  Thus  we  are  motivated  to  work  in  the  phase  domain  in  subsequent 
performance  analyses  in  that  the  distortion  HD(e?a)  is  no  longer  present  and  therefore  not 
of  concern. 

In  this  chapter  we  will  use  both  an  analytical  and  an  experimental  approach  to  the 
prediction  of  the  performance  of  backward  difference  numerical  FM  discrimination  in  the 
presence  of  additive  white  Gaussian  noise  (AWGN).  The  experimental  aspect  consists  of 
computer  simulations  where  finite  length  complex  signal  sequences  consisting  of  the  ideal 
DTFM  signal  and  additive  noise  are  generated  and  processed  for  performance 
measurement  purposes.  By  calculating  a  pair  of  statistics  to  be  described,  an  experimental 
performance  measure  is  established  for  comparison  to  the  analytically  predicted  results.  It 
should  be  noted  that  in  this  research,  the  term  "white"  refers  to  the  fact  that  the  noise 
power  spectral  density  is  constant  over  the  Nyquist  frequency  interval.  In  continuous 
time,  it  is  common  to  refer  to  this  noise  as  being  "bandlimited  white"  noise. 
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The  signal  model  to  be  used  for  both  the  analytical  and  experimental  performance 
predictions  is  depicted  in  Figure  (4-1).  The  model  consists  of  the  ideal  complex  DTFM 
signal  sequence  as  described  by  Equations  (3-1)  through  (3-4),  in  summation  with  a 
complex  Gaussian  noise  sequence,  both  of  which  are  represented  in  Cartesian  form.  The 
discrete-time  noise  sequence  is  the  sampled  bandpass  noise  process 

N(nTs)  =  NtinTs)  +j  •  Nq(nTs)  ,  (4-4) 

where  NfnT)  is  the  Gaussian  distributed  in-phase  component  sequence,  and  Nq(nTs )  is  the 
Gaussian  distributed  quadrature  component  sequence.  The  random  processes  N^nT)  and 
Nq(nT)  are  real-valued  independent  and  identically  distributed  processes  of  zero  mean  and 
variance  a 2.  This  signal-plus-noise  model  is  typical  of  that  used  in  communication 
systems  analyses  [1 1,  36], 


*  K ) 


N,(nTs) 


?K) 


T 

) 


>  XjflT,  ) 


Figure  (4-1).  The  model  of  the  ideal  complex  DTFM  signal  in  the  presence  of  bandlimited 
additive  white  Gaussian  noise. 
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Prior  to  sampling,  we  consider  the  power  spectral  density  of  the  noise,  as 

being  uniformly  distributed  over  the  frequency  interval  -W/2  </<  +W/2  ,  and  of  density 
2N0  Watts/Hz  .  The  corresponding  autocorrelation  sequence,  Rm(T),  found  as  the  inverse 

Fourier  transform  of  the  power  spectral  density,  is 


Rnn(x)  =  2No  W' 


sin(7dfx) 

nWx 


(4-5) 


The  zeros  of  this  function  occur  at  integer  multiples  of  W .  Thus  by  ensuring  that  the 
sample  interval  is  selected  to  correspond  to  the  occurrence  of  these  zeros,  the  samples  of 
the  complex  noise  sequence  will  be  uncorrelated.  Since  these  samples  are  of  a  zero  mean 
Gaussian  process,  they  are  also  independent.  For  analysis  purposes,  we  choose  the 
maximum  required  sample  interval  satisfying  the  Nyquist  criterion,  which  corresponds  to 
the  smallest  non-zero  time  interval  that  yields  uncorrelated  samples.  The  corresponding 
discrete-time  power  spectral  "density"  is  [12] 


Sm(ejQ)  =  2NqFs  ,  -tc  <  Q  <  +7t ,  Q  =  co7’s  .  (4-6) 

Thus  after  sampling,  the  noise  power  is  uniformly  distributed,  and  therefore  considered  to 
be  "white",  over  the  Nyquist  interval.  With  a  spectral  density  of  2N0  Watts/Hz  over  the 
frequency  interval  -W/2  </<  +W/2  Hz,  we  find  that  the  power  of  the  complex  noise 
waveform  is  2N0W=  2NqFs.  This  corresponds  to  twice  the  variance,  o2,  of  both  the 
in-phase  and  quadrature  components.  The  resultant  noise-contaminated  complex 
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DTFM  signal  is  represented  as  the  random  process 


X+(nTs)  =  Xi(nTs)  +j  •  Xq(nTs) ,  (4-7) 


where 


X,(nTs)  =  i{nTs)  +N,(nTs) 


(4-8) 


and 


Xq(nTs)  =  q{nTs)  +Nq(nTs)  .  (4-9) 

Again,  i(nT)  and  q(nTs)  are  as  in  Equations  (3-2)  and  (3-3)  respectively,  with  the  imposed 
modulation  as  in  Equation  (3-4).  Subsequent  analyses  and  simulations  will  be  based  on 
the  established  signal  model  of  Equations  (4-7)  through  (4-9). 

4.1.1  The  Phase  Noise  Probability  Density  Function 

Consider  the  instantaneous  phase  message  estimate 

$/(«r,)=Arctan{^>,  (4-10) 

defined  in  an  analogous  fashion  to  Equations  (3-12)  and  (4-2).  (For  notational  purposes, 
the  tilde  symbol,  ~,  over  a  phase  variable  will  be  indicative  of  noise-contaminated  phase.) 
Of  interest  is  the  resulting  random  phase  noise  sequence,  r\(nTs),  a  constituent  of  the 
random  instantaneous  phase  sequence 
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$,(«n)  =  g[a(nTs )  +  r\(nTs)  +  2n  •  Re(nTs )]  ,  (4-11) 

where  RfnT ’s)  is  some  unknown  sheet  sequence  error,  and  the  true  phase,  a(nTs),  is 

a (tiTs)  =  2nfenTs  +  §(nTs)  -  0  .  (4-12) 

The  phase  noise  sequence  is  defined  relative  to  the  ideal  noise-free  phase  as 

ri(«r,)  =  g[U"Ts)  -  a (nTs)]  .  (4-13) 

Equivalently  the  phase  noise  can  be  defined  relative  to  the  ideal  noise-free  instantaneous 
phase  as 


H(nT.)  =  g[UnTs) -g[a(«r,)]]  .  (4-14) 

Geometrically,  this  process  is  depicted  on  the  phase  plane  in  Figure  (4-2).  Note  that  the 
signal  vector,  s+ ,  can  lie  on  any  sheet,  but  for  simplicity  is  shown  on  sheet  0. 

By  virtue  of  the  independence  of  the  random  variables  X.  and  Xq ,  the  joint 
probability  density  function  (pdf)  of  Xj  and  Xq  is  easily  found  as  the  product  of  the 
corresponding  Gaussian  marginal  pdfs.  After  transforming  the  resulting  joint  pdf  from 
Cartesian  to  polar  coordinates  and  integrating  over  the  envelope  variable,  Middleton  [3, 
37]  has  shown  the  marginal  pdf  of  the  phase  noise  to  be 

/>Cn)  =  h  ■  exp(-p)  +  ^S(T1)  ■  exp(-p  sin2(ri))  •  Erf  (J2p  cos(ri))  .  (4-15) 
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Here,  the  function  Erf  (u)  is  the  error  function,  typically  defined  in  communication  system 
analysis  as 


Erf  (u)  =  j  j=--  exp(-f  )dx  .  (4-16) 

-00  * 


Figure  (4-2).  The  phase  plane  representation  of  the  phase  noise,  T| ,  at  an  arbitrary  time, 
t  =  nTs. 


The  input  (i.e.,  prior  to  phase  estimation)  carrier-to-noise  ratio,  p  ,  is  defined  as 


(4-17) 
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Figure  (4-3)  is  a  plot  of  this  phase  noise  pdf  (herein  referred  to  as  the  Middleton 
distribution)  for  various  values  of  p.  Note  that  as  a  result  of  defining  the  phase  noise 
relative  to  the  true  phase,  a,  we  have  that 

— 7C  <  rj  <  +7t  .  (4-18) 

It  is  also  readily  apparent  that  since  the  input  noise  sequence  is  independent  from  sample 
to  sample,  the  output  (i.e.,  after  phase  estimation)  phase  noise  sequence  will  also  be 
independent  from  sample  to  sample.  Therefore  the  resulting  phase  noise  power  spectral 
density  will  also  be  uniform  over  the  Nyquist  interval,  —  </<  ^  .  This  is  an  additional 
motivation  for  basing  performance  comparisons  in  the  phase  angle  domain  rather  than  the 
phase-difference  domain,  the  latter  of  which  results  in  correlated  consecutive  samples. 

We  observe  that  as  the  carrier-to-noise  ratio  increases,  r|  decreases  such  that  we 
can  use  the  approximations 


sin2(r|)  =  r|2 ,  and  cos(r()  =  1 . 

In  this  case,  the  Middleton  distribution  of  Equation  (4-15)  approaches  a  zero  mean 
Gaussian  distribution  of  variance  (a/a)2.  At  the  other  extreme,  for  values  of  p 
approaching  zero,  the  Middleton  distribution  approaches  that  of  a  uniform  distribution  in 
the  interval  (-7t,7t]. 


85 


1  .  8 


Figure  (4-3).  The  Middleton  probability  density  function  of  the  phase  noise,  r\,  at  various 
values  of  the  carrier-to-noise  ratio,  p. 


4.1.2  The  Phase  Cycle-Slip 

In  addition  to  the  additive  phase  noise,  r|,  which  is  present  in  our  instantaneous 
phase  estimate,  we  are  also  concerned  with  errors  in  the  determination  of  the  phase  plane 
sheet  on  which  the  signal  vector  lies.  In  this  research,  we  will  not  be  concerned  with 
locally  constant  errors  in  the  determination  of  this  integer  value,  rather,  we  will  be 
concerned  with  changes  in  this  error  between  consecutive  sample  intervals.  (Note  that  it 
is  these  changes  in  this  error  sequence  which  lead  to  errors  in  the  resulting  backward 
difference  instantaneous  frequency  sequence.)  A  change  in  the  error  of  the  phase  plane 
sheet  determination  over  a  sample  time  interval  will  be  referred  to  as  a  phase  cycle-slip. 
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Recall  that  our  true  phase  angle  is  as  in  Equation  (4-12), 


a (nTs)  =  2%fenTs  +  fy(nTs )  -  0  , 

corresponding  to  the  phase  of  the  noise-free  signal  vector,  s+(nTs).  Expressed  in  terms  of 
the  phase  plane  sheet  sequence,  r(nT),  and  the  sheet  0  residual  phase  constituents,  we 
have 


a (nTs)  =  g[a(wjTs)]  +  2k  ■  r(nTs )  ,  r  e  integers  .  (4-19) 

As  described  in  Chapter  2,  the  integer  sequence,  r(nTs),  identifies  the  sheet  number  on 
which  the  signal  vector  lies.  Note  that  a  change  in  r  from  time  t=[n-\]Ts  to  time  t  =  nTs 
does  not  constitute  a  phase  cycle-slip.  It  is  simply  representative  of  the  progression  of  the 
signal  vector  across  the  +/-K  boundary  of  the  phase  plane,  during  this  time  interval. 

We  will  be  interested  in  identifying  the  estimate  of  the  sequence  r(nTs),  derived 
from  the  phase  message  estimate,  and  knowledge  of  the  true  phase  sequence,  ct(nTs).  The 
baseline  discriminator  process  applied  to  the  noise-contaminated  DTFM  signal  results  in 
the  instantaneous  frequency  estimate 


VfrTd  =  F'.gEfa(flT')-%&n-\\T.yL ,  (4-20) 

with  (j);(/?7y)  defined  as  in  Equation  (4-10).  In  likeness  to  Equation  (4-3)  we  can  define 
our  noise-contaminated  phase  message  estimate  as 
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(4-21) 


3U«rs)  =  Ts  •  { i  %'d{kTs) + m-t,)}  , 

k=  0 

A 

which  will  be  used  to  identify  the  corresponding  sheet  sequence  estimate,  R(nTs).  Thus  a 
phase  cycle-slip  occurs  whenever  there  is  a  change,  from  time  t  =  [n-X\Ts  to  t  =  nTs,  in  the 
integer  valued  error  sequence, 

ReinTs)  =  R(nTs)  -  r(nTs )  .  (4-22) 

4.1.3  The  Phase  Noise  and  Phase  Cycle-Slip  Performance  Statistics 

To  predict  the  performance  of  the  discriminator  phase  message  estimate  of 
Equation  (4-21)  via  simulation,  it  is  necessary  to  introduce  measurement-derived  statistics 
which  are  representative  of  the  phase  noise  and  phase  cycle-slip  errors.  This  section 
presents  the  simulation  technique  used  to  determine  the  phase  noise  sample  sequence  from 
knowledge  of  the  true  phase  sequence.  Subsequently,  the  sheet  sequence  estimate  is 
presented,  which  can  be  used  to  determine  phase  cycle-slip  occurrences. 

Observe  that  in  the  simulation  process,  the  true  phase  as  described  by  Equation 
(4-12)  is  known,  since  it  is  needed  to  generate  the  in-phase  and  quadrature  pairs, 
{i(nTs),q(nTs)},  for  the  finite  set  of  indices  n  =  0,1,2,...,(NPTS-1).  For  simplicity,  each 
vector  component  generated  by  the  sine  and  cosine  functions  is  left  as  normalized  to  unity, 
i.e.,  a  =  1 .  Each  simulation  run  begins  with  the  addition  of  this  finite-length  complex 
signal  sequence  and  a  newly  generated  complex  Gaussian  noise  sequence  of  zero  mean 
and  variance  2o2.  The  new  noise  sequence,  {N^nT^N^nT)},  n  =  0,1,2,...,(NPTS-1), 
generated  for  each  simulation  run,  has  the  properties  described  in  Section  4. 1,  namely  : 
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1 )  Nf  and  Nq  are  zero  mean,  Gaussian  uncorrelated  random  variables  of  variance  o2  at  any 
particular  time,  t  =  nTs,  and 

2)  {NfnT),Nq{nT^}  and  {Nl(mT),Nq(mT)}  are  uncorrelated  for  all  indices  with  the 
exception  of  m  =  n. 

The  addition  of  the  signal  and  noise  finite  length  sequences  results  in  the 
noise-contaminated  sequence,  {Xfn Ts),Xq(n Ts) } ,  n  =  0,1,2,...,(NPTS-1),  generated  at  the 
beginning  of  a  simulation  run,  having  a  fixed  carrier-to-noise  ratio  equal  to  l/(2a2). 
During  a  simulation  run,  the  real  sequence 


d(nTs)  =  Ts  •  %'d(nTs) ,  n  =  0, 1, 2, ...,  (NPTS-  1)  (4-23) 

is  generated  using  a  backward  difference  numerical  FM  discrimination  method.  This 
allows  us  to  calculate  our  phase  message  estimate  simply  as 

$d(nTs)=i  d(kTs),  n  =  0,1,2,..., (NPTS -1).  (4-24) 

k=0 

Consider  first,  the  determination  of  a  phase  noise  performance  measure.  Just  as 
the  zero-mean  Gaussian  distribution  is  completely  characterized  by  its  variance,  the 
Middleton  phase  noise  pdf  in  Equation  (4-14)  is  completely  characterized  by  its  variance, 
F{r|}.  The  phase  noise  variance  is  a  decreasing  function  of  p,  and  is  bounded  at  the  upper 
extreme  by  the  variance  of  a  uniform  distribution  in  (-%, n],  and  at  the  lower  extreme  by 
the  variance,  (p/a)2,  of  a  Gaussian  distribution.  Thus  we  are  interested  in  an  estimate  of 
the  variance  of  the  phase  noise,  as  one  of  the  performance  measures. 
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Knowledge  of  the  true  phase,  a(nTs),  allows  us  to  consider  determining  a  sample 
function  of  the  phase  noise  sequence 

r\(nTs)  =  g&t(nTl)-g{a(nTs)Ji 

-g[^d(nTs)-a(nTs)]  .  (4-25) 

Here,  we  have  taken  advantage  of  the  fact  that 

U’>T,)=dtA>T,)]  (4-26) 


as  shown  in  Appendix  B,  and  the  modulo  2n  arithmetic  property  presented  in  Equation 
(3-49).  From  a  practical  viewpoint,  this  means  that  we  can  determine  the  finite  length 
phase  noise  sequence,  r| (nT),  from  Equation  (4-25)  as 

r\(nT,)  =  gEihinTs)  -  a (tiTs)] ,  n  =  0, 1, 2, ...,  (NPTS- 1)  .  (4-27) 

As  an  estimate  of  the  variance  of  the  phase  noise,  we  define  the  performance  statistic 


(4-28) 


(4-29) 
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With  the  ability  to  obtain  a  phase  noise  sample  function  and  knowledge  of  the  true  phase, 
we  now  consider  the  phase  cycle-slip  performance  measurement. 

In  Appendix  B,  we  have  represented  the  noise-contaminated  phase  message 
estimate  as 


$d(nTs)  =  a (nTs)  +  r](nTs)  +  2tc  •  Re(nTs )  ,  (4-30) 

where  Re{nTs )  is  as  described  in  the  appendix,  an  integer  valued  sheet  sequence  error. 

From  Equation  (4-30)  and  the  phase  noise  sequence  of  Equation  (4-27),  knowledge  of  the 
true  phase  sequence  allows  us  to  identify  the  sheet  sequence  estimate  as 

2:r  •  R(nTs)  =  {$rf(«7V)  -  -n(»n)}  -gElhinTs)  -  ri(«7;)] ,  n  =  0, 1, 2, ...,  (NPTS-  1)  . 

(4-31) 


In  practice,  the  calculation  in  Equation  (4-31)  is  found  to  be  sensitive  to  computer 
representation/roundoff  errors  when  a  is  an  odd  multiple  of  +/-k,  and  can  lead  to  +/-2iz 
errors  when  used  in  conjunction  with  Equation  (4-22)  to  determine  a  sample  function  of 
the  sheet  sequence  error,  Re(nTs).  Therefore,  although  Equations  (4-22)  and  (4-3 1)  yield 
further  insight  into  the  nature  of  the  phase  sheet  sequence  error,  we  select  an  alternate 
method  of  determining  Re(nTs).  An  equivalent  calculation  for  practical  purposes  is 

2ti  •  Re(nTs)  =  $rf(A?r,)  -  a{nTs)  - y\(nTs) ,  n  =  0, 1, 2, ...,  (NPTS -  1)  ,  (4-32) 
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directly  from  Equation  (4-30).  By  counting  the  number  of  negative  transitions,  NCOUNT, 
and  positive  transitions,  PCOUNT,  in  the  error  sequence  found  in  Equation  (4-32),  we 
have  a  measure  of  phase  cycle-slip  performance  for  each  simulation  run.  Alternatively,  we 
can  form  the  ratio 


a  _  NCOUNT  +  PCOUNT  TCOUNT  (4-33) 

Pcs  ~  NPTS  MTS  ’  V  * 

as  an  estimate  of  the  probability  of  a  phase  cycle-slip  occurring  at  any  particular  sample  of 
the  phase  message  estimate.  Here,  TCOUNT,  is  the  total  number  of  transitions,  either 
positive  or  negative. 

Taken  together,  the  phase  cycle-slip  count  and  the  phase  noise  variance  estimate 
form  a  set  of  measures  which  can  be  used  to  predict  the  performance  of  the  baseline 
numerical  FM  discrimination  process.  This  technique  is  not  unlike  that  of  Rice's  [6],  in 
which  both  "cycle-skips"  and  output  signal-to-noise  ratio  analytical  measures  are  used  to 
summarize  continuous-time  FM  discriminator  performance.  However,  in  a  similar  fashion 
to  Develet's  analysis  of  angle  demodulation  [20,  39],  we  have  chosen  to  stay  in  the 
phase-domain  with  both  performance  measures. 
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4.2  Recovery  of  the  Sheet  Sequence,  r{nTs),  in  the  Presence  of  AWGN 

We  have  previously  introduced  methods  of  obtaining  phase  noise  variance  and 
phase  cycle-slip  performance  measures  for  the  backward  difference  numerical  FM 
discrimination  technique,  via  experimental  computer  simulations.  In  these  simulations, 
noise  sample  function  sequences  obtained  from  a  Gaussian  random  process  are  added  to 
the  ideal  complex  DTFM  signal,  for  subsequent  analyses.  We  will  now  be  concerned  with 
the  theoretical  phase  cycle-slip  performance  prediction  in  more  detail,  for  comparison  to 
the  simulation  results.  The  noise-contaminated  signal  model  is  as  presented  in  Section  4. 1 . 

Consider  the  accumulation  process  of  Equation  (4-21),  which  is  used  as  our 
noise-contaminated  phase  message  estimate.  From  this  we  can  write 

h(nTs)  =  -  l]r,)  +  Ts  •  tyd{nTs) 

=  $,([»  -  l]r,)  +g[$/(«r,)  -  $,([«  -  1]7V>] •  (4-34) 

Using  Equation  (4-26)  and  modulo  2k  arithmetic,  Equation  (4-34)  is  rewritten  as 

%d(nTs)  =  $«,([»  -  1]^)  +g\g[fa(nTs )]  -  1]^)]] 

=  $«/([»-  l]T.)+g!&d(nT,)-$4ln-l]T,)\  .  (4-35) 

After  application  of  Equation  (4-30)  and  modulo  2k  arithmetic  in  Equation  (4-35), 
rearrangement  of  terms  yields 

r\(nTs)  =  v(nTs)  -  2k  •  { Re(nTs ) -Re([n  -  1]^)} ,  (4-36) 
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where 


v(nTs)  =  r|([«  -  l]Ts)  +g[a'd(nTs )  +  r| d(nTs)\  -  a'd(nTs )  (4-37) 

The  terms  a d(nTs)  and  r\d(nTs)  are  defined  as  the  backward  differences 

a'd(nTs)  =  <x(nTs) - a([n  -  l]r,)  (4-38) 


and 


t|  'd(nTs)  =  r\  (nTs)  -  ri([»  -  1J77,).  (4-39) 

Equation  (4-36)  is  a  key  result  in  that  it  allows  us  to  determine  the  probability  of  a  phase 
cycle-slip  occurring  at  a  particular  time  t  =  nTs.  This  equation  indicates  that  when  v(nTs) 
exceeds  n  in  magnitude,  a  phase  cycle-slip  identified  as 

2ti  •  Rcs(nTs)  =  2;r  •  { Re(nTs )  -  Re{[n  -  1J77,)}  (4-40) 

will  occur  as  a  "compensation"  such  that  equality  with  r|(«rs)  holds.  The  result  is  that  the 
phase  noise,  r\(nTs),  remains  on  sheet  0  of  the  modulation  referenced  phase  plane  as 
defined,  and  as  indicated  in  Equation  (4-18). 

Figure  (4-4)  is  representative  of  the  various  pdfs  which  arise  at  a  particular  time, 
t  =  nTs,  in  the  formation  of  the  random  process  v(nT)  .  For  simplicity  we  use  the  notation 

a’d=  a'pTX  and  =  T]([w  - 1]^),  -  r\(nTs),  v  =  v(nT),  for  these  random  variables. 
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Figures  (a)  through  (c)  are  for  the  case  where  the  quantity  (r|_  -  a'J  is  greater  than  zero. 
Similarly,  figures  (d)  through  (f)  are  for  the  case  where  the  quantity  (r|  -  a'd)  is  less  than 
zero. 


(d) 


(e) 


(0 


Figure  (4-4).  Representative  pdfs  arising  in  the  formation  of  v(«7J. 


Case  1 :  (r)  -  a'J  >  0 

a)  The  pdf  of  Cn  -  r|  +  a’d)  given  r\  and  a’d . 

b)  The  pdf  of  g[r|  -  r\ _  +  a'J  given  r|  and  a’d . 

c)  The  pdf  of  v  =  r|-  a’d  +  g[rj  -  T]+  a'J  given  r)  and  a’d ,  pjvlrj  ,a'J  . 
Case  2  :  (r|  -  a'J  <  0 

d)  The  pdf  of  (n  -  r|_  +  a'J  given  r\  and  a’d. 

e)  The  pdf  of  g[ri  -  r\  +  a'J  given  r\  and  a’d . 

f)  The  pdf  of  v  =  ti  -  a’d  +  g[q  -  ri_  +  a'J  given  ri_  and  a'd ,  /?<(v|r|_  ,a'J  . 
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From  the  representative  pdfs  of  figures  (c)  and  (f),  it  is  possible  to  determine  the 
probability  of  a  phase  cycle-slip  under  various  modulation  conditions.  The  three  cases  of 
interest  will  be  the  unmodulated  carrier  with  a  perfectly  tuned  receiver,  the  unmodulated 
carrier  with  an  off-tuned  receiver,  and  finally,  the  angle  modulated  carrier  with  an 
off-tuned  receiver.  The  following  sections  will  cover  each  of  these  cases,  presenting  both 
the  theoretically  predicted  and  the  corresponding  simulation  results. 

4.2.1  Case  I :  The  Unmodulated  Carrier  and  Perfectly  Tuned  Receiver 

We  begin  with  the  case  where  the  carrier  is  unmodulated  and  there  exists  no  tuning 
error,  i.e.,  the  case  where  cn(nTs)  =  0.  In  this  case,  a 'd  =  -9  ~(-6)  —  0.  Consider  the 
determination  of  the  probability  that  v  exceeds  7t,  which  corresponds  to  the  occurrence  of 
a  positive  phase  cycle-slip  compensation,  2uRcs,  at  some  particular  time  t  —  nTs.  We 
represent  this  probability  as  ppcs  =  Pr{v  >  tc}.  Similarly  we  represent  the  probability  of  a 
negative  phase  cycle-slip  as  pncs  —  Pr{v  <  -7t}.  The  probability  of  a  phase  cycle-slip 
occurrence  (either  positive  or  negative)  is  simply 


Pcs  ~Pncs  +  Ppcs  ■  (4  41) 

This  is  true  regardless  of  the  modulation  scenario,  since  the  occurrence  of  a  positive  phase 
cycle-slip,  the  occurrence  of  a  negative  phase  cycle-clip,  and  the  occurrence  of  no  phase 
cycle-slip,  form  a  set  of  mutually  exclusive  and  exhaustive  events. 

Due  to  the  symmetry  of  the  Middleton  phase  noise  pdf,  and  due  to  the  fact  that  the 
backward  difference,  a'd ,  is  equal  to  zero,  it  is  readily  apparent  that  the  probabilities  of  a 
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negative  phase  cycle-slip  and  of  a  positive  phase  cycle-slip  are  the  same.  In  Figure  (4-4  c), 
the  shaded  area  represents  the  probability  of  a  positive  phase  cycle-slip,  given  a  specific 
value,  r|_  ,  of  the  random  variable  r|  .  Integrating  over  all  possible  values  of  r|  we  obtain 


U  JC+T 1— 


Pncs  =  Ppcs  =  Pr{ v  >  7C>  =J  J  p  >(vlri_)  p(r\-)dv  dx\- 


(4-42) 


Note  that  integration  on  the  variable  r\  is  from  0  to  n  rather  than  -k  to  n,  since  for  values 
of  t)  less  than  zero,  Pr{v>7i}  =  0.  To  simplify  Equation  (4-42)  we  employ  the  fact  that 

1=g[v]  .  (4-43) 

as  implied  by  Equation  (4-36).  Thus  values  of  v  in  the  range  n  to  (7t  +  r\_ )  map  to  values 
of  r|  in  the  range  -7i  to  (-71  +  r\  ).  Taking  advantage  of  this  mapping,  we  rewrite  Equation 
(4-42)  as 


7t+T|- 


Pncs  =  Ppcs  =  Pr{-7I  <  T|  <  -71  +  T1_ }  =  J /?(t)_)  J  p(r\)dr\  dr \- 


(4-44) 


Note  that  both  x\  and  r|  are  Middleton  distributed  as  in  Equation  (4-15).  Given  the  input 
carrier-to-noise  ratio,  p  ,  we  can  use  Equation  (4-44)  to  determine  the  probability  of  a 
positive  phase  cycle-slip,  which  in  this  special  case  is  the  same  as  the  probability  of  a 
negative  phase  cycle-slip. 
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Before  presenting  a  comparison  of  the  theoretical  and  the  experimental 
performance  predictions  for  various  input  carrier-to-noise  ratios,  we  consider  a  specific 
example  in  which  the  carrier-to-noise  ratio  is  allowed  to  approach  zero.  In  this  example, 
the  pdfs  of  r|_  and  rj  approach  that  of  a  uniform  distribution  in  the  interval  (-7t,7t].  From 
Equation  (4-44)  we  obtain  the  probability  of  a  negative  phase  cycle-slip  as 


I 

8 


We  conclude  that  with  perfectly  tuned  reception  of  an  unmodulated  carrier  in  the  presence 
of  AWGN,  the  backward  difference  numerical  FM  discriminator  will  result  in  a  negative 
phase  cycle-slip  occurring  at  a  particular  time,  t  =  nTs,  with  a  probability  approaching  1/8 
as  p  approaches  zero.  As  previously  indicated,  the  symmetry  of  this  example  implies  that 
the  probability  of  a  positive  phase  cycle-slip  also  approaches  1/8.  From  Equation  (4-41), 
the  probability  of  a  phase  cycle-slip  therefore  approaches  1/4. 

Figures  (4-5  a)  through  (4-5  d)  present  both  theoretical  and  experimental  results 
serving  as  prediction  of  the  performance  of  the  backward  difference  discriminator,  with 
regard  to  phase  cycle-slip  occurrences.  Throughout  the  various  simulations,  the  sequence 
length,  NPTS,  is  equal  to  1024  samples. 
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Figure  (a)  presents  the  analytically  predicted  number  of  negative  phase  cycle-slips 
per  NPTS  samples,  as  p  is  varied  from  -20  to  +40  dB  in  2  dB  increments.  This  predicted 
count  is  found  as  the  product  of  pncs  and  NPTS.  Here,  pnC5  is  found  by  numerical 
integration  of  the  integral  of  Equation  (4-44).  Figure  (b)  presents  the  corresponding 
experimental  results,  which  are  arrived  at  using  the  technique  outlined  in  Section  4.1.3. 
Rather  than  presenting  the  results  in  the  form  of  a  predicted  probability  as  in  Equation 
(4-33),  we  have  simply  plotted  the  variable  NCOUNT  versus  the  carrier-to-noise  ratio. 
Note  that  the  experimental  results  agree  quite  well  with  the  theoretical  results.  It  is  also 
noted  that  as  p  decreases,  the  experimental  negative  phase  cycle-slip  count  approaches 
values  near  1024/8  =  128,  as  expected. 

Similarly,  figures  (c)  and  (d)  respectively  present  the  theoretically  and 
experimentally  predicted  number  of  positive  phase  cycle-slips  in  a  sequence  of  length  1024 
samples.  The  simulations  presented  in  figure  (d)  plot  the  variable  PCOUNT  versus  the 
input  carrier-to-noise  ratio.  As  expected,  the  performance  predictions  are  the  same  as  for 
the  negative  phase-cycle  slip  measures,  in  the  currently  treated  case  of  the  unmodulated 
carrier  and  perfectly  tuned  receiver. 

Our  attention  now  turns  to  the  case  where  there  exists  a  tuning  error,  fe,  in  our 
DTFM  signal.  This  error  is  associated  with  the  receive  end  of  our  communication  system, 
which  is  modeled  as  the  quadrature  downconverter  (or  its  equivalent)  and  the  backward 
difference  discriminator.  The  error  mechanism  is  therefore  referred  to  as  an  off-tuned 
receiver. 
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Figure  (4-5  a).  Analytically  predicted  number  of  negative  phase  cycle-slips  vs.  p 
Case  I :  The  unmodulated  carrier  and  perfectly  tuned  receiver. 


Figure  (4-5  b).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  I :  The  unmodulated  carrier  and  perfectly  tuned  receiver. 
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Figure  (4-5  c).  Analytically  predicted  number  of  positive  phase  cycle-slips  vs.  p  ; 
Case  I :  The  unmodulated  carrier  and  perfectly  tuned  receiver. 


Figure  (4-5  d).  Experimentally  predicted  number  of  positive  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  I :  The  unmodulated  carrier  and  perfectly  tuned  receiver. 
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4.2.2  Case  II :  The  Unmodulated  Carrier  and  Off-tuned  Receiver 

We  now  analyze  the  case  where  the  carrier  is  unmodulated,  but  there  exists  a 
tuning  error  of  fe  Hz.  In  this  case,  a’d  =  (2%fenTs  -  0)  -  (27 ife  [n  -  Y\TS  -  0)  =  2nfeTs. 
Consider  the  determination  of  the  probability  that  v  is  less  than  -n,  which  corresponds  to 
the  occurrence  of  a  negative  phase  cycle-slip  compensation,  2nRcs,  at  some  particular  time, 
t  -  nTs.  As  before,  we  represent  this  probability  as  pncs  =  Pr{v  <  -7t}.  This  corresponds 
to  the  scenario  where  (tj  -  <x'd)  is  less  than  zero,  as  shown  in  Figure  (4-4  f).  Integrating 
over  all  contributing  values  of  r\_  we  obtain 


For  the  outer  integral  of  Equation  (4-45),  values  of  r|  in  the  range  [-n,a'd]  will  result  in 
non-zero  contributions  over  the  inner  integral.  Thus  the  limits  of  integration  for  the  outer 
integral  are  established. 

Likewise,  for  the  scenario  where  (r|  -  a'd)  is  greater  than  zero  corresponding  to  a 
positive  phase  cycle-slip,  we  write 

Ppcs  =J  j  /?>(vlrt_) -p(r\-)dv  dr\-  (4-46) 

* 

For  the  outer  integral  of  Equation  (4-46),  only  values  of  t|  in  the  range  [cl'^k]  will 
contribute  to  the  overall  integration. 
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To  simplify  the  integrals  of  Equation  (4-45)  and  (4-46)  we  again  employ  Equation 
(4-43).  For  (q_  -  a'd)  less  than  zero,  values  of  v  in  the  range  -n  to  -(7C+q_-a'd)  map  to 
values  of  q  in  the  range  (rc+ri  -oc'd)  to  n.  This  allows  us  to  rewrite  Equation  (4-45)  as 


j  P(n)dr\  dr |_  . 

("+n--a^) 


(4-47) 


Similarly,  for  (q_  -  a’d)  greater  than  zero,  values  of  v  in  the  range  %  to  (K+r\  -a'd)  map  to 
values  of  q  in  the  range  —k  to  (tc^  q  -cx.'d).  Equation  (4-46)  is  rewritten  as 

It  -(n-T|--HXrf) 

Ppcs  =  J  p(y\-)  J  p{r\)dr\  dr\-  .  (4_48) 

ad 

As  in  Case  I,  we  determine  the  probability  of  a  phase  cycle-slip  (either  positive  or 
negative)  using  Equation  (4-41).  It  is  interesting  to  note  that  by  employing  the  symmetry 
of  the  Middleton  pdf,  a  simple  change  of  variable  in  Equation  (4-47)  along  with  a  change 
in  the  order  of  integrations  yields  the  property 


Pncs(CLd,  p)  —Ppcsir&d,  P)  • 


(4-49) 


In  the  above  expression,  we  have  explicitly  indicated  the  dependence  of  these  probabilities 
on  both  the  tuning  offset,  a'd  =  2nfeTs ,  and  the  input  carrier-to-noise  ratio,  p. 


Before  presenting  a  comparison  of  the  Case  II  theoretical  and  experimental 
performance  predictions  for  various  input  carrier-to-noise  ratios,  we  consider  a  specific 
example  in  which  the  carrier-to-noise  ratio  is  allowed  to  approach  zero.  As  stated 
previously,  the  pdfs  of  ti_  and  n  approach  that  of  a  uniform  distribution  in  the  interval 

In  this  example,  we  also  allow  the  tuning  offset  to  approach  the  positive  extreme 
value,  a'd  =  -be.  From  Equation  (4-47)  we  obtain  the  probability  of  a  negative  phase 
cycle- slip  under  this  condition  as 


2k  2% 


7b 

-~~x  f  (7C  -  ri_)z/n- 
4k2  4 


We  conclude  that  with  maximally  positive  off-tuned  reception  of  an  unmodulated  carrier  in 
the  presence  of  AWGN,  the  backward  difference  numerical  FM  discriminator  will  result  in 
a  negative  phase  cycle-slip  occurring  at  a  particular  time,  /  =  nTs,  with  a  probability 
approaching  1/2  as  p  approaches  zero.  At  the  same  time,  from  Equation  (4-48)  the 
probability  of  a  positive  phase  cycle-slip  is  zero.  From  the  property  presented  in  Equation 
(4_49)5  this  example  implies  that  when  a'd  =  -k,  the  probability  of  a  positive  phase 
cycle-slip  also  approaches  1/2  as  p  approaches  zero.  Likewise,  from  Equation  (4-47)  the 
probability  of  a  negative  phase  cycle-slip  is  zero,  when  a'd  =  -k. 
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Figures  (4-6  a)  through  (4-6  d)  present  both  theoretical  and  experimental  results 
serving  as  prediction  of  the  performance  of  the  backward  difference  discriminator,  with 
regard  to  phase  cycle-slip  occurrences  for  Case  II.  Throughout  the  various  simulations, 
the  sequence  length,  NPTS,  is  again  equal  to  1024  samples.  The  tuning  offset  in  these 
figures  is  a'd  =  2nfeTs  =  7t/10. 

Figure  (a)  presents  the  analytically  predicted  number  of  negative  phase  cycle-slips 
per  NPTS  samples,  as  p  is  varied  from  -20  to  +40  dB  in  2  dB  increments.  This  predicted 
count  is  found  as  the  product  of  pncs  and  NPTS.  Here,  pncs  is  found  by  numerical 
integration  of  the  integral  of  Equation  (4-47).  Figure  (b)  presents  the  corresponding 
experimental  results,  which  are  arrived  at  using  the  technique  outlined  in  Section  4.1.3. 

As  in  Case  I,  rather  than  presenting  the  results  in  the  form  of  a  predicted  probability  using 
Equation  (4-33),  we  have  simply  plotted  the  variable  NCOUNT  versus  the  carrier-to-noise 
ratio.  Note  that  the  experimental  results  agree  quite  well  with  the  theoretical  results. 

Similarly,  figures  (c)  and  (d)  respectively  present  the  theoretically  and 
experimentally  predicted  number  of  positive  phase  cycle-slips  in  a  sequence  of  length  1024 
samples.  The  simulations  presented  in  figure  (d)  plot  the  variable  PCOUNT  versus  the 
input  carrier-to-noise  ratio.  As  expected,  the  presence  of  a  positive  tuning  offset  increases 
the  chance  of  a  negative  phase  cycle-slip,  and  decreases  the  chance  of  a  positive  phase 
cycle-slip. 

Next,  we  will  be  interested  in  the  case  where  there  exists  a  tuning  error,  fe,  in 
addition  to  a  message  signal,  <J)(0  ,  in  our  DTFM  signal.  This  is  the  final  and  most  general 
modulation  case  consistent  with  our  established  DTFM  model,  which  will  be  considered. 
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Figure  (4-6  b).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 


(4  simulation  runs);  Case  II :  The  unmodulated  carrier  and  off-tuned  receiver 

(2tu/7s  =  =  71/1  °)- 
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Figure  (4-6  c).  Analytically  predicted  number  of positive  phase  cycle-slips  vs.  p 
Case  II :  The  unmodulated  carrier  and  off-tuned  receiver  (2nf  T  =  =  n/10). 
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Figure  (4-6  d).  Experimentally  predicted  number  of positive  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  II :  The  unmodulated  carrier  and  off-tuned  receiver 
(2t ifeTs  =  %5e  =  7i/10). 


4.2.3  Case  lU  :  The  Angle  Modulated  Carrier  and  Off-tuned  Receiver 

The  final  case  to  be  considered  is  the  general  case  in  which  there  exists  both  a 
tuning  error  and  angle  modulation  imposed  on  the  carrier.  Demodulation  of  the  noise-free 
DTFM  signal  in  this  case  yields  o.'JnTs)  =  2nfeTs  +  §(nTs)  -  <)>([«  -  1]^).  We  will  build 
upon  the  results  of  Equations  (4-47)  and  (4-48)  to  arrive  at  phase  cycle-slip  probabilities 
for  the  present  modulation  case. 

It  is  observed  that  in  the  noise-contaminated  DTFM  model  which  we  have 
adopted,  the  distribution  of  phase-difference  amplitude  values,  a'd,  is  statistically 
independent  of  the  distribution  of  each  of  the  phase  noise  random  variables,  r|  and  rj_. 

This  amplitude  distribution  will  be  represented  as  Pj(y),  where  the  variable  Y  signifies  any 
of  the  possible  values  of  the  quantity  y  =  a’d  at  any  arbitrary  time,  t  =  nTs.  As  previously 
established,  these  amplitude  values  range  from  -7t  to  +7t.  Given  the  amplitude 
distribution,  pfy>),  the  probability  of  a  negative  phase  cycle-slip  is  readily  found  as 

Pncs  =  J  pr(y)  j  /?0l-)  I  p(y\)dx\  dr |_  dy  .  (4-50) 

(7t+n--a^) 

Equation  (4-50)  can  be  viewed  as  a  weighted  average  of  Case  II  negative  phase  cycle-slip 
probabilities.  Previous  expressions  for  the  probability  of  a  negative  phase  cycle-slip  are 
simply  special  cases  of  Equation  (4-50),  in  which  the  quantity  >>  =  <x’d  takes  on  a  single 
known  value. 

Similarly  the  probability  of  a  positive  phase  cycle-slip  is 
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(4-51) 


it  it  -(»-n-+Orf) 

Jp(ri-)  J  p(r\)dr\dr\-dy  . 

~n  y  -n 

This  can  likewise  be  interpreted  as  the  weighted  average  of  Case  II  positive  phase 
cycle-slip  probabilities.  Once  again,  previous  expressions  for  the  probability  of  a  positive 
phase  cycle-slip  are  simply  special  cases  of  Equation  (4-51),  in  which  the  quantity  y  =  a'd 
takes  on  a  single  known  value. 

We  now  present  the  specific  example  in  which  the  imposed  angle  modulation  is 
sinusoidal  with 


m{nTs)  =  sin(27 ifmnTs  +  0m), 


(4-52) 


such  that 


<K«r*)  -  <K[w  -  1]^)  =  Ts  ■  hD(nTs)  *  {2kA /•  m(nTs)} .  (4-53) 

In  the  presence  of  a  tuning  offset,  employing  results  from  Appendix  A  yields 

y(nTs)  =  \iy  +aysm(2nfmnTs  +  (4-54) 


where 


[ly  —  2%feTs  , 


(4-55) 


and 
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(4-56) 


2  A/"  *  r  fm  \ 

ay  =  --sm(K--)  . 

Random  sampling  of  the  amplitude  of  y(nTs )  yields  the  distribution 

Py(y) =  7==  >  ~ay  -(y~  M-y)  -  ay  •  (4-57) 

To  carry  out  the  numerical  integration  of  Equation  (4-50)  and  (4-51),  we  perform  a 
weighted  average  of  the  inner  double  integrals,  over  a  finite  set  of  subintervals  spanning 
the  amplitude  interval,  [(-a  +  py),(a  +  p^)].  The  weights  used  are  the  corresponding 
probabilities  of  y  being  in  a  particular  subinterval.  Corresponding  to  each  subinterval,  any 
value  of  the  amplitude,  y,  in  this  subinterval  can  be  used  in  the  calculation  of  the  inner 
double  integrals.  Consider  first  the  numerical  calculation  of  Equation  (4-51),  the 
probability  of  the  occurrence  of  a  negative  phase  cycle-slip. 

The  approximation  to  Equation  (4-50)  which  is  used  to  find  the  probability  of  a 
negative  phase  cycle-slip  is 

2b-\  y  n 

Pncs  =  S  Pk-\ /?Cn-)  j  p(y)dr\  dx\- ,  yk^y  ^+1 ,  (4-58) 

fc=1  (jt+n_-a  j) 

where  the  probability  weights,  P#  are  found  as 

n+i 

Pk  =  J  Pr(y)dy  .  (4-59) 

yk 
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Here,  the  amplitude  subintervals  are  a  result  of  quantizing  the  modulating  signal  into  2b 
equally  spaced  discrete  values  of  y,  where  b  is  some  positive  integer.  (In  the  simulations 
to  be  presented,  b  =  5.)  The  subintervals  are  \yk ,  y^x  ],  k  =  1,2,3, ...,(26-l),  with  yx 
corresponding  to  (~a  +  py)  and  y2h_1  corresponding  to  (a  +  p^).  In  the  sinusoidal 
modulation  example,  with  pjy)  defined  as  in  Equation  (4-57),  the  weights  become 

Pk  =  i ■  [sin-1  (  )  -  sm  ■(  ^  ),  (4-60) 

For  a  numerically  calculated  upper  bound  on  Equation  (4-50)  we  employ  Equation  (4-58) 
with  y=ykH  ,  since  Pr{v  <  -7t}  is  an  increasing  function  of  y.  Similarly,  a  lower  bound 
on  the  probability  of  a  negative  phase  cycle-slip  is  determined  using  Equation  (4-58),  by 
setting  y=yk . 

Likewise,  an  approximation  to  the  probability  of  the  occurrence  of  a  positive  phase 
cycle-slip  described  by  Equation  (4-51)  is 

2 »_!  7i  -<*— n-+<4) 

Ppcsz'Z  Pk  \p(r 1-)  J  p(v)dr]dr\-  ,  yk<y<yk+i  .  (4-61) 

k=  1  y  -n 

The  weights,  Pk,  remain  as  in  Equation  (4-59).  For  a  numerically  calculated  upper  bound 
on  Equation  (4-51),  we  employ  Equation  (4-61)  with  y=yk,  since  Pr{v  >  7i}  is  a 
decreasing  function  of  y.  Similarly,  a  lower  bound  on  the  probability  of  a  positive  phase 
cycle-slip  is  determined  using  Equation  (4-61)  by  setting  y  =yk+l  . 
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Figures  (4-7)  and  (4-8)  present  both  the  theoretically  and  experimentally  predicted 
phase  cycle-slip  performance  measures,  for  two  distinct  Case  III  modulation  scenarios. 

The  first  scenario  is  presented  in  Figures  (4-7  a)  through  (4-7  d),  with  the  relative  message 
sampling  rate,  y,  set  at  17.067,  and  the  DTFM  message  index,  8m,  set  at  0.2930.  Figure 
(a)  presents  the  upper  and  lower  bounds  on  the  analytically  predicted  number  of  negative 
phase  cycle-slips.  The  corresponding  experimental  results  are  presented  in  figure  (b),  and 
are  in  agreement  with  the  analytically  predicted  performance.  Similarly,  figure  (c)  presents 
the  analytically  predicted  positive  phase  cycle-slip  occurrences,  while  figure  (d)  presents 
the  experimental  results.  Again,  the  analytical  and  experimental  performance  predictions 
are  in  agreement  with  each  other. 

Likewise,  Figures  (4-8  a)  through  (4-8d)  present  the  theoretical  and  experimental 
phase  cycle-slip  performance  predictions  for  the  modulation  scenario  where  y  remains  set 
at  17.067,  but  8m  is  changed  to  0.5859.  As  expected,  increasing  the  DTFM  modulation 
index  has  the  effect  of  increasing  the  probability  of  a  phase  cycle-slip.  Figure  (a)  presents 
the  upper  and  lower  bounds  on  the  analytically  predicted  number  of  negative  phase 
cycle-slips.  The  corresponding  experimental  results  are  presented  in  figure  (b),  and  are  in 
agreement  with  the  analytically  predicted  performance.  Similarly,  figure  (c)  presents  the 
analytically  predicted  positive  phase  cycle-slip  occurrences,  while  figure  (d)  presents  the 
experimental  results.  As  seen  in  these  figures,  the  analytical  and  experimental  performance 
predictions  are  in  agreement  with  each  other. 


112 


Figure  (4-8  a).  Analytically  predicted  number  of  negative  phase  cycle-slips  vs.  p  ; 

Case  III :  The  angle  modulated  carrier  and  off-tuned  receiver  (2%feTs  =  n8e  =  n/10,  with 
sinusoidal  modulation,  y  =  17.067,  8m  =  0.5859). 


Figure  (4-8  b).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  III :  The  angle  modulated  carrier  and  off-tuned  receiver 
( '2nfeTs  =  n8e  =  k/10,  with  sinusoidal  modulation,  y  =  17.067,  5m  =  0.5859). 


4.3  The  Mean-Square  Phase  Noise  Performance  Measure,  EP 

The  emphasis  in  this  chapter  thus  far  has  been  on  the  phase  cycle-slip  performance 

measure.  Our  attention  now  turns  toward  the  measured  mean-square  phase  noise,  as  a 
phase  noise  variance  estimate,  for  performance  prediction  purposes.  Recall  from 
Equations  (4-27)  through  (4-28)  that  this  measure  requires  the  knowledge  of  the  true 
phase  sequence,  a (nT),  as  in  the  phase  cycle-slip  measure.  Taken  together,  both 
measures  are  used  to  predict  the  performance  of  the  backward  difference  numerical  FM 
discriminator  in  the  presence  of  additive  white  Gaussian  noise. 

We  begin  by  considering  the  special  case  in  which  the  imposed  angle  modulation  is 
sinusoidal  as  in  Equation  (4-52).  In  the  presence  of  a  tuning  error,  fe,  noise-free  phase 
recovery  yields 

a (nTs)  =  Z  y(kTs)  =  2nfenTs  +  p  sm(2n/mnTs  +  -  f)  -  0  ,  (4-62) 

where  P  is  A/ /  fm  as  defined  in  Equation  (3-73).  Note  that  the  "information-bearing" 
phase  signal  component  is  the  term 

<$>(nTs)  =  p  sin(2Tif mnTs  +  Qm  -  f) .  (4-63) 

The  output  noise-contaminated  phase  is  from  Equation  (4-30)  represented  as 

$d(«7V)  =  2KfenTs  +  §(nTs)  -  0  +  q (nTs)  +  2n  ■  Re(nTs ) ,  (4-64) 
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which  motivates  us  to  consider  an  output  signal-to-noise  ratio  performance  measure  using 
our  knowledge  of  §(nTs),  and  the  phase  noise  variance  estimate. 

The  output  phase  signal-to-noise  ratio  (for  sinusoidal  modulation)  is  defined  as 


1  v  =  _fe_ 

V{r] }  7  2V{r\) 


(4-65) 


The  factor  y  in  the  above  equation  represents  the  effect  of  an  ideal  lowpass  filter  of 
bandwidth  fm  Hz,  applied  to  the  recovered  phase  signal.  Thus  the  effect  of  the  output 
phase  noise,  t],  is  measured  over  the  message  bandwidth.  The  factor  p2/2  is 
representative  of  the  power  of  the  sinusoidal  modulating  signal. 

Figure  (4-9)  presents  the  results  of  both  the  analytically  and  experimentally 
predicted  output  phase  SNR  measure,  as  a  function  of  the  input  carrier-to-noise  ratio,  p. 
Figure  (a)  presents  the  analytical  results,  while  figure  (b)  presents  the  results  of  a  single 
simulation  run.  For  the  experimental  performance  prediction,  the  numerically  found  phase 
noise  variance,  F{r|},  has  been  replaced  by  its  estimate,  Ep ,  using  a  simulated  DTFM 
sequence  of  length  NPTS  =  1024  samples.  In  both  figures,  the  modulation  index,  8m, 
takes  on  the  values  0.0586,  0. 1172,  0.2930,  and  0.5860,  with  y  fixed  at  17.067.  (These 
values  of  8m  correspond  to  P  values  of  1,  2,  5  and  10  respectively,  for  the  current  value 
of  y.)  Note  that  the  experimentally  predicted  measure  is  in  close  agreement  with  that  of 
the  analytically  predicted  measure,  for  the  given  modulation  indices.  It  is  interesting  to 
observe  that  a  single  curve  from  Figure  (4-9)  is  sufficient  to  serve  as  phase  noise 
performance  prediction;  the  parameter  P  merely  shifts  (increases)  the  ordinate  by  the 
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Figure  (4-9  a).  Analytically  predicted  output  phase  SNR  in  the  message  bandwidth; 
(sinusoidal  modulation,  y  =  17.067). 
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Figure  (4-9  b).  Experimentally  predicted  output  phase  SNR  in  the  message  bandwidth; 
(sinusoidal  modulation,  y  =  17.067). 


amount  101og{P2/2}.  Note  also,  however,  that  as  the  phase  noise  pdf  approaches  a 
uniform  distribution  (for  values  of  p  less  than  0  dB  or  so)  the  output  phase  SNR  measure 
loses  its  meaning.  In  this  limiting  situation,  the  phase  noise  very  effectively  randomizes 
the  recovered  instantaneous  phase  and  thus  eliminates  the  signal.  Due  to  the  fact  that  the 
phase  noise  at  time  t  =  nTs  is  modeled  as  being  statistically  independent  of  the  phase  noise 
at  time  t  =  [»- 1  ]  Ts,  the  backward  difference  does  not  help  to  alleviate  this  effect.  (In  fact 
as  indicated  in  [6],  Rice  [38],  Stumpers  [5],  and  Middleton  [4],  have  shown  that  for 
the  continuous-time  case,  the  message  signal  is  reduced  by  the  factor  (1  -  e'p ).) 

We  make  the  observation  that  for  a  fixed  sample  rate,  Fs,  and  fixed  message 
bandwidth,  fm,  the  frequency  deviation,  A/,  can  be  increased  such  that  for  a  given  input 
carrier-to-noise  ratio,  the  output  phase  signal-to-noise  ratio  can  be  enhanced.  This  is 
done,  however,  at  the  expense  of  increasing  the  probability  of  phase  cycle-slip 
occurrences,  as  seen  in  Figure  (4-10).  In  comparison  to  Figures  (4-7)  and  (4-8),  we  see 
that  increasing  the  DTFM  message  index,  8m,  to  0.8203  (by  increasing  A f)  gives  rise  to 
additional  phase  cycle-slips,  at  input  carrier-to-noise  ratios  as  high  as  20  dB. 

A  more  realistic  scenario  is  that  the  sample  rate  is  decreased  in  order  to  facilitate 
processing  on  a  computing  device.  For  a  fixed  message  bandwidth  and  deviation 
frequency,  the  result  is  that  y  decreases  while  8m  increases,  both  towards  unity.  Again  this 
is  done  at  the  expense  of  additional  phase  cycle-slips. 
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Figure  (4-10  a).  Analytically  predicted  number  of  negative  phase  cycle-slips  vs.  p  ; 
Case  III :  The  angle  modulated  carrier  and  ofT-tuned  receiver  (2nfeTs  =  tz5  =  tc/10,  with 
sinusoidal  modulation,  y  =  17.067,  5m  =  0.8203). 


Figure  (4-10  b).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  III :  The  angle  modulated  carrier  and  off-tuned  receiver 
(2%feTs  -  n5  =  7t/10,  with  sinusoidal  modulation,  y  =  17.067,  8m  =  0.8203). 


Figure  (4-10  c).  Analytically  predicted  number  of  positive  phase  cycle-slips  vs.  p  ; 
Case  III :  The  angle  modulated  carrier  and  off-tuned  receiver  {2%feTs  =  k8  =  71/ 10,  with 
sinusoidal  modulation,  y  =  17.067,  5m  =  0.8203). 


Figure  (4-10  d).  Experimentally  predicted  number  of  positive  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Case  III :  The  angle  modulated  carrier  and  off-tuned  receiver 
(2%feTs  -  =  7t/10,  with  sinusoidal  modulation,  y  =  17.067,  8m  =  0.8203). 
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4.4  Relating  Instantaneous  Frequency  Aliasing  and  the  Phase  Cycle-Slip 

In  the  absence  of  noise,  we  have  previously  established  in  Equation  (3-67)  the 
requirement  to  prevent  instantaneous  frequency  aliasing  due  to  overmodulation.  With 
a'piT)  defined  as  in  Equation  (4-38),  this  requirement  is  restated  as 

—7t  <  a d(nTs)  <  +7T  .  (4-66) 

With  noise  as  modeled  in  the  present  research,  the  new  requirement  becomes 

—7i  <  ud(nTs)  +  r\d(nTs )  <  +7t .  (4-67) 

Note  that  even  in  the  latter  noise-contaminated  case,  the  requirement  of  Equation  (4-66) 
remains.  Thus  it  is  best  to  refer  to  violation  of  Equation  (4-66)  as  overmodulation. 
Violation  of  Equation  (4-67)  will  be  referred  to  as  instantaneous  frequency  aliasing. 

Further  consideration  of  the  phase  cycle-slip  as  derived  from  Equation  (4-37)  leads 
us  to  conclude  that  analysis  of  the  term 

g[ad(nTs)  +  r\d(nTs)\  , 

will  uniquely  determine  when  cycle-slips  occur.  More  specifically,  a  phase  cycle-slip  is 
avoided  under  the  condition 

ga'/pT.) + n>r,)] = a >r.) + n  >r,)  (4-68) 
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This  condition,  however,  is  simply  a  restatement  of  Equation  (4-67).  Thus  the  phase 
cycle-slip  probabilities  previously  determined  have  the  alternate  interpretation  of  being  the 
probabilities  of  instantaneous  frequency  aliasing  events.  Aliasing  across  the  +7t  boundary 
is  the  same  as  a  negative  phase  cycle-slip,  while  aliasing  across  the  -71  boundary  is  the 
same  as  a  positive  phase  cycle-slip. 

4.5  Concluding  Remarks 

As  already  indicated,  by  taking  into  account  both  the  occurrences  of  phase 
cycle-slips  and  the  output  phase  SNR  measure,  we  have  a  method  of  predicting  the 
performance  of  the  backward  difference  discriminator.  In  the  previous  section  we  have 
not  addressed  the  fact  that  the  occurrence  of  a  phase  cycle-slip  is  dependent  on  the  phase 
noise  sequence  sample  function,  such  that  these  distinct  noise  contributions  are  not 
independent  of  each  other.  In  this  regard,  these  two  noise  effects  should  not  be 
considered  directly  additive.  Rather,  careful  consideration  must  be  used  to  properly 
combine  the  effects,  if  so  desired.  Currently,  we  will  consider  it  to  be  sufficient  to  leave 
the  performance  measures  as  separate  quantities. 
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CHAPTER  5 


Simulation  of  a  Method  of  Discriminator  Performance  Enhancement 
for  Efficient  Processing  of  Large  P  DTFM  Signals  at  Low  Input  CNR 


5.1  Objective 

We  have  established  in  Chapter  4  the  representation  of  the  recovered  phase  from 
the  backward  difference  discriminator,  in  units  of  radians,  as 

$</( nTs)  =  a (nTs)  +  r](nTs)  +  2nRe(nTs) ,  (5-1) 

where 

a (nTs)  =  2nfetiTs  +  §(nTs)  -  0  .  (5-2) 

Here,  2 itfenTs  is  a  ramp  function  resulting  from  the  tuning  error ,/„  <\>(nT)  is  the  phase 
message  signal,  0  is  an  arbitrary  but  constant  reference  phase,  and  r\(n7s)  is  the  Middleton 
distributed  phase  noise  process.  The  term  2nRe(nTs)  is  representative  of  the  occurrence  of 
phase  cycle-slips  and  has  the  properties  outlined  in  Chapter  4  and  Appendix  B.  This 
recovered  phase  representation  is  the  result  of  modeling  the  modulated  input  signal  as  the 
Cartesian  sum  of  our  complex  DTFM  signal  with  complex  Gaussian  distributed  noise  of 
uniform  spectral  density  over  the  Nyquist  frequency  band.  In  this  research  focus  will  be 
on  the  reduction  of  phase  cycle-slip  occurrences,  by  reducing  the  modulation  index  of  the 
noise-contaminated  DTFM  signal. 

Of  interest  to  us  is  the  modulation  scenario  in  which  the  original  continuous-time, 
constant  envelope,  frequency  modulated  signal  is  characterized  by  a  "large"  Af  /fm  ratio 
(i.e.,  large  P).  Currently  we  will  consider  the  term  "large"  to  imply  that  the  minimum 
required  sample  rate,  after  conversion  of  the  input  to  the  generalized  pre-envelope  form,  is 
such  that  the  message  signal  amplitude  changes  by  no  more  than  10  to  20%  between 
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consecutive  samples.  (Recall  that  the  message  signal,  m(nTs),  is  normalized  in  amplitude.) 
Note  that  for  efficient  and/or  real-time  processing  of  such  a  DTFM  signal,  it  is  desirable  to 
perform  computations  at  as  low  a  sample  rate  as  possible.  However,  as  we  have  seen  in 
previous  results,  lowering  the  sample  rate  increases  the  probability  of  the  occurrence  of 
phase  cycle-slips  for  the  given  resultant  value  of  relative  message  sample  rate,  y.  In  effect, 
by  lowering  the  sample  rate,  the  DTFM  modulation  index,  5,  increases  towards  its 
maximum  allowable  value  of  unity,  beyond  which  the  signal  is  considered  to  be 
over-modulated. 

As  indicated  in  Figures  (4-10  a)  through  (4-10  d),  we  find  that  the  phase  cycle-slip 
probability  increases  substantially  with  increases  in  modulation  index.  This  can  be  seen  in 
comparing  these  results  to  those  of  Figures  (4-8  a)  through  (4-8  d).  At  a  10%  tuning 
offset,  by  increasing  8m  from  .5859  to  .8203  the  phase  cycle-slip  count  per  1024  samples 
increases  from  approximately  zero,  to  tens  of  counts  at  values  of  input  CNR  in  the  range 
10  to  20  dB.  In  this  chapter,  we  will  develop  and  evaluate  a  method  of  enhancing  the 
performance  of  the  backward  difference  numerical  FM  discriminator,  using  frequency 
feedback  techniques.  In  particular,  by  simulating  the  noise-contaminated  complex  DTFM 
signal  as  in  Chapter  4,  observed  phase  cycle-slip  counts  and  phase  noise  variance 
estimation  will  be  used  to  compare  the  enhanced  demodulator  to  the  established  backward 
difference  discriminator  baseline. 

Special  note  is  made  of  the  fact  that  in  this  chapter,  for  convenience,  instantaneous 
frequency  estimates  are  left  in  units  of  normalized  radian  frequency,  rather  than  scaling  by 
Fs  to  convert  to  units  of  radians  per  second. 
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5.1.1  The  FM  with  Feedback  (FMFB)  Demodulator 

The  continuous-time  FMFB  technique  as  proposed  by  Chaffee  [IS]  and  furthered 
by  Enloe  [19]  and  others,  was  originally  introduced  as  a  means  of  reducing  the  FM 
threshold  effect  for  large  deviation,  low  carrier-to-noise  ratio  FM  systems.  The  basic 
block  diagram  representation  of  the  FMFB  device  is  presented  in  Figure  (5-1).  The 
related  and  better  known  phase-lock  loop  (PLL)  device  is  shown  by  Develet  [20]  to  be 
equivalent  to  the  FMFB  demodulator  under  the  correct  conditions. 


Bandpass 
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- > 
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Filter 


(LPF) 


Figure  (5-1).  The  continuous-time  (analog)  FM  with  Feedback  (FMFB)  demodulator. 


As  explained  by  Enloe,  both  the  PLL  and  the  FMFB  demodulators  use  the  a  priori 
information  that  the  ratio  of  frequency  deviation  to  message  signal  bandwidth,  A/ /fm,  is 
large.  Referring  to  Figure  (5-1),  the  function  of  the  Voltage  Controlled  Oscillator  (VCO) 
is  to  track  the  slowly  varying  message-signal-induced  frequency  excursions.  At  the  same 
time,  the  Bandpass  Filter  (BPF)  responds  only  to  a  much  narrower  band  of  noise  centered 
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around  these  frequency  excursions.  The  function  of  the  Analog  Limiter/Discriminator  and 
Lowpass  Filter  (LPF)  is  to  provide  the  necessary  information  to  the  VCO  to  continuously 
track  these  excursions.  Note  that  an  error  in  the  VCO  frequency  relative  to  the  input 
waveform  will  result  in  a  proportionate  control  voltage,  acting  in  a  negative  feedback 
manner,  such  that  the  error  is  reduced  and  continuously  maintained  near  zero.  Thus  the 
discriminator  is  presented  a  narrower  band  of  noise  so  that  any  resulting  phase  noise  is 
substantially  lowered  at  the  discriminator  output. 

5.1.2  The  Numerical  FM  with  Feedback  Demodulator 

We  arrive  at  a  form  of  the  numerical  FMFB  demodulator  by  converting  to 
discrete-time  the  analog  FMFB  device  of  Figure  (5-1).  By  additionally  incorporating  the 
generalized  pre-envelope  representation,  we  arrive  at  the  system  shown  in  Figure  (5-2). 
The  input  to  the  device  is  the  noise-contaminated  generalized  pre-envelope  sequence, 
X+(nTs,fo).  The  data  registers  shown  are  simply  indicative  of  the  fact  that  all  necessary 
calculations  for  each  output  estimate  must  be  performed  within  the  sample  clock  interval, 
Ts  seconds.  It  is  assumed  that  the  instruction  clock  of  the  programmable  device  used  to 
implement  the  process  and  the  sample  clock  are  derived  from  a  common  clock  source  such 
that  an  integer  number  of  instruction  cycles  are  available  between  consecutive  time 
samples. 
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Figure  (5-2).  A  numerical  implementation  of  the  FMFB  demodulator. 


For  a  description  of  its  operation,  we  start  at  time  t  -  nTs  at  the  input  to  the  device. 
The  function  of  the  Complex  Multiply  module  is  to  reduce  the  modulation  index  of  the 
input  sequence  by  reducing  the  frequency  deviation,  A/.  With  the  input  sequence 
represented  as 

X+(nTs,fo)  =  U+(«r,)l  •  exp{/'[27t/0«rs  +(J)(«7Ti)-0  +t|(h7’j)]},  (5-3) 

the  resultant  bandpass  filtered  product  is 

{ X+(nTs,f0 )  •  exp{-jhP(nTsJ0)}  }bpf  =  \A+(nTs)BPF\  •  exp{j$e(nTs)} ,  (5-4) 

or  expressed  in  Cartesian  form, 


\A+(tiTs)bpf\  •exp{/$e(w7,J)}  =  {Xel(nTs)  +JXeq(nTs)}  =Xe(nTs).  (5-5) 

The  notation  {.}BPF  is  representative  of  the  effect  of  the  complex  bandpass  filter  in  the 
device,  shown  as  the  BPF  module.  (This  filter  processes  the  result  of  the  Complex 
Multiply  module.)  Here,  \A+(rtT^\  is  the  noise-perturbed  envelope  sequence  associated 
withX+(«7^s,/j.  Likewise,  \A+{nT)BPF  \  is  the  resulting  envelope  after  bandpass  filtering. 
With  the  definition 


a+(nTs)  =  K-l  {hlp(kTs)  *y([k- l]r,)} ,  (5-6) 


the  phase  prediction  sequence  is 
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bP(nTs,f0)  =  2; ifonT,  +  &  +(nTs)  =  2%  [fa -fe]nTs  +  a +(nTs), 


(5-7) 


where  hlp(nTs )  is  the  impulse  response  of  the  phase  predictor  lowpass  filter.  We  have  also 
used  the  fact  that 


27 ifo  =  27 l  \fo  -/«]. 


(5-8) 


The  resulting  residue  modulation  becomes 

%e{nTs)=g[{a(nTs)  +  r[(nTs)-a+(nTs)}BPF\  .  (5-9) 

The  intent  of  the  bandpass  filtering  operation  is  to  pass  a  narrow  band  of  frequency 
excursions  centered  around  zero  hertz,  and  to  reject  the  noise  energy  outside  of  this  band. 
At  the  same  time,  this  filter  must  remain  wide  enough  in  frequency  such  that  the  residue 
modulation  accurately  represents  the  reduced  index  phase  modulation.  Employing  the 
backward  difference  discriminator,  the  Numerical  FM  Discriminator  module  processes  this 
reduced  modulation  index  signal,  producing  an  estimate  of  the  original  modulating  signal. 
Although  not  shown,  this  estimate  can  subsequently  be  processed  further,  using  for 
example  lowpass  filtering,  to  enhance  the  output  signal-to-noise  ratio.  Thus  available  to 
the  Phase  Predictor  module  for  generating  the  current  phase  prediction,  are  all 
instantaneous  frequency  estimates  prior  to  the  current  time,  t  =  nTs.  With  the  above 
calculations  complete,  a  new  sample  interval  begins  with  a  sample  clock  transition  to 
t  =  [tj+1]T,  and  the  entire  process  continues. 
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The  shown  form  of  the  phase  predictor  consists  of  a  fixed  feedback  gain  factor,  K, 
the  lowpass  filter,  LPF,  and  an  accumulator  which  acts  as  an  integrator  inverting  the 
differentiation  process  of  the  numerical  discriminator.  Injection  of  the  center  frequency 
estimate  after  scaling  by  the  sample  interval,  effectively  sets  the  nominal  operating 
frequency  of  the  Numerically  Controlled  Oscillator  (NCO)  module.  The  shown  NCO 
allows  for  quantization  of  the  feedback  phase  prediction  such  that  the  required  cosine  and 
sine  trigonometric  functions  can  be  readily  implemented  using  look-up  table  techniques. 
The  number  of  bits  used  in  the  quantizer  determines  both  the  size  of  the  look-up  table,  and 
the  amplitude  resolution  of  the  instantaneous  phase  prediction.  The  number  of  bits  can  be 
chosen  to  produce  a  resolution  which  is  commensurate  with  the  maximum  change  between 
consecutive  samples  expected  in  the  modulation  sequence  . 

For  practical  reasons,  the  output  of  the  accumulator  in  the  Phase  Predictor  module 
is  retained  in  a  modulo  2k  fashion,  bounding  the  accumulated  result  in  the  interval  (-71,71] 
and  easing  quantizer  implementation.  Overall,  with  proper  choice  of  feedback  gain, 
complex  bandpass  filtering  and  phase  prediction  lowpass  filtering,  the  numerical 
implementation  acts  in  an  analogous  manner  to  its  continuous-time  counterpart  and 
effectively  tracks  instantaneous  frequency  excursions. 
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5.2  The  Reconstituted  Numerical  FM  with  Feedback  (RNFMFB)  Demodulator 

We  now  introduce  a  form  of  demodulator  device  to  be  referred  to  as  the  Type  I 
RNFMFB  demodulator,  shown  in  Figure  (5-3).  This  device  operates  just  as  the  previously 
described  numerical  FMFB  demodulator  of  Figure  (5-2).  Flowever,  from  this  previous 
description  it  is  observed  that  we  have  available  to  us  the  instantaneous  phase  prediction 
sequence, 


g[bP(nTs,f0)\  =  g[2nf0nTs  +  a,{nTs)}  ,  (5-10) 

which  has  been  used  to  reduce  the  modulation  index  of  the  original  input  sequence, 
X+(nTs,f0).  Since  the  message  output  is  no  longer  taken  directly  from  the  numerical 
discriminator,  we  modify  Equation  (5-6)  as 

d +(nTs)  =  K-t  { hlp(kTs )  *  g[de([k-  1]T,)]} ,  (5-1 1) 


where 


de{nTs )  =  $e(pTs)  -  $,([«  -  1]T,)  . 


(5-12) 


It  is  apparent  that  the  rate  of  change  of  the  instantaneous  phase  prediction  sequence  of 
Equation  (5-10)  represents  the  frequency  compression  that  has  been  effected  upon  the 
original  modulated  signal.  Thus  with  the  appropriate  time  alignment  as  shown  in  the 
Phase-Difference  Reconstitution  module,  this  compression  signal  can  be  recombined  with 
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Figure  (5-3).  The  Type  I  reconstituted  numerical  FM  with  feedback  demodulator. 


the  recovered  residue  modulation  estimate,  reconstituting  the  original  modulation.  The 
time  alignment  takes  into  account  the  group  delay  of  the  complex  bandpass  filter  function, 
BPF.  The  output  is  now  identified  as  this  reconstituted  instantaneous  frequency  estimate. 

The  Type  I  RNFMFB  device,  interesting  in  its  own  right,  is  left  for  future  research 
efforts.  Currently  we  emphasize  the  analysis  of  a  specific  implementation  of  the  Type  II 
RNFMFB  demodulator,  as  presented  in  Figure  (5-4).  The  Type  II  device  is  distinguished 
from  the  Type  I  in  that  the  signal  fed  back  to  the  Phase  Predictor  module  is  now  the 
reconstituted  message  estimate,  rather  than  the  residue  constituent,  g[ de(nTs)\.  Note  that 
for  the  implemented  RNFMFB  demodulator,  as  shown  in  the  figure,  both  the  bandpass 
and  lowpass  filter  functions  have  been  omitted  for  computational  efficiency  and  reduction 
of  complexity.  By  additionally  designing  the  acquisition  and  demodulation  process  such 
that  Ts  is  as  large  as  permissible  for  proper  message  recovery,  we  increase  the  number  of 
instructions  available  between  consecutive  samples. 

In  practice,  we  find  that  cost  effective  programmable  Digital  Signal  Processing 
(DSP)  devices  are  characterized  by  their  ability  to  perform  the  operations  of 
multiply-accumulate  and  table  look-ups  quite  efficiently  in  terms  of  executed  instructions. 
However,  these  devices  typically  cannot  perform  calculations  such  as  division  and 
arctangent  in  an  efficient  manner.  As  a  consequence  of  these  device  characteristics,  the 
presented  RNFMFB  demodulator  becomes  attractive  to  the  communications  engineer, 
given  that  the  extra  processing  requirements  can  be  implemented  in  a  small  percentage  of 
the  total  set  of  instructions  required  for  demodulation.  We  therefore  potentially  gain  more 
instructions  by  increasing  Ts  than  are  invested  in  the  additional  processing.  In  the  absence 
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Numerical  FM  Discriminator 
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Figure  (5-4).  Implemented  Type  II  reconstituted  numerical  FM  with  feedback  demodulator. 


of  this  processing  advantage,  analysis  of  the  performance  of  the  RNFMFB  device  relative 
to  the  established  baseline  enhances  basic  understanding  of  the  operation  of  numerically 

implemented  frequency  feedback  techniques. 

To  facilitate  the  description  of  the  operation  of  the  RNFMFB  demodulator,  we 

refer  to  the  phase-domain  model. 

5.2.1  The  Implemented  RNFMFB  Demodulator  Phase-Domain  Model 

The  phase-domain  model  of  the  implemented  RNFMFB  demodulator  is  shown  in 
Figure  (5-5),  where  we  have  made  use  of  the  z-transform  time  delay  notation,  z  . 

Referring  to  the  figure,  we  find  that 

g[bP(nTsjo)]  =  g[2itf0nTs  +  0.5  iy([k-  1]7Y)]  (5-13) 

such  that 

$,2(nTs)  =  g[d2(nTs,f0)Jo] 

=  g[g[fyp(nTs,f0)]  -g[<M[w  -  1  ]Ts,fo)]  ~  2tc/o] 

=s[0.5K[»-iF,)1 .  <5-14> 

Since 

y(nTs)  =  $[(nTs)  +  %'2(nTs) ,  (5-15) 
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we  find  that  the  message  output  has  the  property 


-2k  <y(nTs )  <  +2k  ,  (5-16) 

which  allows  us  to  simplify  Equation  (5-14)  and  write 

$/2(«r,)  =  0.5j>([«-  l]r,)  .  (5-17) 


Figure  (5-5).  The  phase-domain  model  of  the  implemented  Type  II  RNFMFB 
Demodulator. 
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With  the  output  of  the  secondary  discriminator  established,  we  are  now  interested  in 
determining  an  expression  for  the  primary  discriminator  output.  First  we  note  that 


$e(nTs)  =  g\g[2%f0nTs  +  a  (nTs)  +  Ti(«r,)]  -g[kp(nTsJ0)\] 

=  g[a(/i7’,)  +  ri(w7,,)-0.5  Z  y([k- IF,)] ,  (5-18) 

where  we  have  made  use  of  Equations  (5-2)  and  (5-13).  From  the  above  and  Equation 
(5-8)  we  can  write 


^(nTs)=g[a/d(nTs)+r]'AnTs)-0.5y([n-l]Ts)]  ,  (5-19) 

where  ad(nTs)  and  r\d(nTs)  are  as  defined  in  Equations  (4-38)  and  (4-39).  Equivalently, 
(nTs)  =  a'd(nTs)  +  x\'d(nTs)  -  0.5y([«  -  IF)  +  2: iR*cs(nTs),  (5-20) 

where  the  integer  sequence,  R*s(nTs),  is  required  to  ensure  that  we  retain  the  property 

-n  <  $i (nTs)  ^  •  (5-21) 

This  integer  sequence  is  representative  of  the  occurrence  of  phase  cycle-skips  in  the 
primary  discriminator  of  the  phase-domain  model,  and  has  properties  that  are  unique  to 
the  RNFMFB  device.  Thus  when  the  quantity 

a'd(nTs)  +  y\'d(nTs)-0.5y({n-\]Ts) 
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exceeds  7t  in  magnitude,  a  phase  cycle-skip  occurs  as  compensation  to  ensure  validity  of 
Equation  (5-21).  Combining  the  primary  and  secondary  discriminator  outputs  as  in 
Equation  (5-15),  we  find  from  Equations  (5-17)  and  (5-20)  that  the  message  output 
becomes 


fcnTs)  =  a'd(nTs)  +  x\'d(nTs)  +  2t iR*cs(nTs)  .  (5-22) 

Note  that  Equation  (5-22)  implies  that  the  output  phase  noise  is  identical  to  that  which 
would  have  been  obtained  using  the  backward  difference  discriminator  alone.  Therefore, 
in  subsequent  simulations,  we  expect  the  predicted  output  phase  noise  variance 
performance  measure  to  be  consistent  with  the  baseline  discriminator.  At  best  then,  we 
only  expect  enhanced  performance  via  reduced  phase  cycle-skip  occurrences. 

Observe  that  with  regard  to  the  implemented  RNFMFB  demodulator,  the  intent  of 
the  device  is  to  divide  the  discrimination  process  equally  between  two  numerical  FM 
demodulators  identified  as  the  primary  and  secondary  discriminators.  A  consequence  of 
the  fact  that  we  are  processing  large  P  DTFM  signals,  is  that  the  message  signal  changes 
in  amplitude  by  only  a  small  percentage  between  consecutive  samples.  This  in  turn  results 
in  a  large  relative  message  sampling  rate,  y.  This  a  priori  information  allows  us  to  use  the 
current  instantaneous  frequency  estimate  as  the  predicted  value  of  the  next  estimate.  With 
a  feedback  factor  of  K  =  0.5,  the  frequency  deviation  of  the  signal  presented  to  the 
primary  discriminator  is  potentially  reduced  by  50%.  We  also  note  that  with  this  value  of 
K,  phase  cycle-skips  are  confined  to  the  primary  discriminator. 
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In  the  absence  of  a  previous  phase  cycle-skip  at  time  t  -  [n-\~\Ts,  the  quantity 


controlling  the  occurrence  of  a  cycle-skip  becomes 

a>r„)  +  rfo/iT.)  -  0  5 Ml" -  IF,)  =  a'/nT.)  -  0.5 o!d([n  -  IF,)  +  e'X'F),  (5-23) 

where 

e 'd(nTs)  =  n >T.)  -  0.5r!' ([//  -  1]^)  .  (5-24) 

At  values  of  input  carrier-to-noise  ratio,  p  =  a2/(2o2),  greater  than  10  dB  or  so,  the  phase 
noise  is  approximately  Gaussian,  with  a  variance  of  less  than  (o/a)2.  For  P  >  4,  this  leads 
to  pre-filtered  output  phase  signal-to-noise  ratios  greater  than  20  dB.  This  in  turn  implies 
that  the  primary  discriminator  noise  term  identified  in  Equation  (5-24)  is  small  relative  to 
the  reduced  index  signal, 


a[(nTs)  =  a'd(nTs)  -  0.5a'd([n -  IF,).  (5'25) 

Using  the  pessimistic  modulation  scenario  in  which  the  message  signal  is 
sinusoidal,  the  relative  message  sampling  rate  can  be  used  to  bound  the  message  prediction 
error.  It  can  be  shown  that 


A/M  max 


=  2  cos[tc  •  ( 


y-i 

2y 


)] 


(5-26) 


represents  this  maximum  message  amplitude  change  between  consecutive  samples. 
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5.3  RNFMFB  Demodulation  Performance  Simulation  Results 


The  Type  II  RNFMFB  demodulator  of  Figure  (5-4)  has  been  implemented  for 
comparison  to  the  performance  of  the  established  backward  difference  numerical 
discriminator  baseline.  Each  demodulator  is  presented  four  (for  a  total  of  eight) 
independently  generated  sample  functions  of  the  noise-contaminated  DTFM  process,  of 
length  1024  samples.  Modulation  parameters  have  been  selected  which  demonstrate  the 
typical  performance  of  these  demodulators  for  large  (3  systems.  Specifically,  the  relative 
message  sample  rate  is  set  at  y  =  17.067  at  a  10%  tuning  offset,  with  a  message  index  of 
5m  =  .8203.  This  choice  of  parameters  coincides  with  a  P  value  of  14. 

These  parameters  were  chosen  to  ensure  that  the  maximum  message  amplitude 
change  between  consecutive  samples  as  indicated  in  Equation  (5-26),  is  less  than  20%. 
This  is  consistent  with  our  clarification  of  the  term  "large  P".  In  the  presented  simulations, 
the  selected  y  value  yields  a  maximum  message  amplitude  change  of  18.4%.  The  selected 
parameters  also  represent  a  progression  from  the  performance  results  previously  given  in 
Chapter  4  for  the  backward  difference  discriminator. 

Once  again,  we  present  both  the  probability  of  occurrence  of  phase  cycle-slips  and 
phase  noise  variance  estimates  as  performance  measures.  However,  special  note  is  made 
with  regard  to  the  phase  cycle-slip  counts.  As  already  indicated,  the  RNFMFB 
demodulator  produces  a  phase  cycle-skip  sequence  which  can  potentially  change  by  more 
than  2tc  between  consecutive  samples.  This  fact  has  been  taken  into  account  in  the 
simulations  to  be  presented,  such  that  phase  cycle  changes  greater  than  2 n  are  specifically 
looked  for  during  the  counting  process.  In  the  simulations  presented,  no  such  phase  cycle 
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changes  were  observed  at  or  above  0  dB  input  carrier-to-noise  ratio.  Therefore,  we 
properly  refer  to  the  phase  cycle  changes  observed  for  the  RNFMFB  demodulator,  as 
phase  cycle-slips. 

Figures  (5-6  a)  through  (5-6  d)  present  the  observed  number  of  both  negative  and 
positive  phase  cycle-slips  for  the  backward  difference  and  RNFMFB  demodulators.  In 
particular  we  note  the  significantly  enhanced  performance  of  the  RNFMFB  device 
indicated  in  Figure  (5-6  b)  over  that  of  the  baseline  discriminator  results  in  Figure  (5-6  a). 
For  values  of  input  carrier-to-noise  ratio  between  10  and  20  dB,  the  RNFMFB  device 
produces  no  observed  cycle-slips,  whereas  the  non-enhanced  backward  difference 
demodulator  produces  tens  of  cycle-slips. 

Estimates  of  phase  noise  variance  were  also  made,  and  converted  to  output  phase 
signal-to-noise  ratios  in  the  message  signal  bandwidth,  as  in  Chapter  4.  In  comparing  the 
results  of  the  backward  difference  discriminator  of  Figure  (5-7  a)  with  the  RNFMFB 
demodulator  of  Figure  (5-7  b),  we  note  no  observable  difference.  This  is  consistent  with 
the  resultant  phase  noise  indicated  in  Equation  (5-22)  for  the  RNFMFB  demodulator,  and 
theoretical  results  obtained  for  the  backward  difference  discriminator  in  Appendix  B. 

We  conclude  that  the  RNFMFB  demodulator  is  a  viable  method  of  performance 
enhancement  for  the  backward  difference  numerical  FM  discriminator,  under  the 
appropriate  modulation  conditions.  Specifically  we  have  demonstrated  the  ability  of  the 
RNFMFB  demodulator  to  operate  at  combined  low  sample  rate  and  low  input 
carrier-to-noise  ratios  for  large  P  systems. 
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Input  Carrier-to-Noise  Ratio  (dB) 

Figure  (5-6  a).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Backward  Difference  Discriminator.  (2nfeTs  =  nd  =  nf  10,  with 
sinusoidal  modulation,  y  =  17.067,  8m  =  0.8203). 


Figure  (5-6  b).  Experimentally  predicted  number  of  negative  phase  cycle-slips  vs.  p 
(4  simulation  runs);  RNFMFB  Demodulator.  (2%feTs  =  k5  =  n/10,  with  sinusoidal 
modulation,  y  =  17.067,  8  =  0.8203). 


Input  Carrier-to-Noise  Ratio  (dB) 

Figure  (5-6  c).  Experimentally  predicted  number  of  positive  phase  cycle-slips  vs.  p 
(4  simulation  runs);  Backward  Difference  Discriminator,  (2%feTs  =  nd  =  7r/10,  with 
sinusoidal  modulation,  y  =  17.067,  5m  =  0.8203). 


Input  Carrier-to-Noise  Ratio  (dB) 

Figure  (5-6  d).  Experimentally  predicted  number  of  positive  phase  cycle-slips  vs.  p 
(4  simulation  runs);  RNFMFB  Demodulator,  (2t ifeTs  =  7tS  =  tc/10,  with  sinusoidal 
modulation,  y  =  17.067,  8m  =  0.8203). 


Figure  (5-7  b).  Experimentally  predicted  output  phase  SNR  in  message  bandwidth; 
RNFMFB  Demodulator  (sinusoidal  modulation,  y  =  17.067,  5m  =  .8203). 
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CHAPTER  6 

Conclusions  and  Areas  of  Future  Research 
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6.1  Conclusions 


In  this  research,  a  performance  analysis  has  been  performed  regarding  a  class  of 
numerical  FM  demodulators  operating  in  the  presence  of  additive,  spectrally  uniform, 
Gaussian  noise.  Both  theoretical  and  computer  simulation  techniques  have  been 
employed,  and  have  given  consistent  performance  prediction  results.  The  class  of 
numerical  demodulators  which  have  been  analyzed  apply  the  backward  difference  estimate 
of  the  time  derivative  of  the  phase,  to  estimate  the  instantaneous  frequency  of  an 
angle-modulated  complex  exponential  sequence.  Several  computationally  different 
approaches  to  backward  difference  numerical  FM  demodulation  have  been  shown  to 
produce  equivalent  results.  Basic  concepts  regarding  the  representation  and  processing  of 
complex-valued,  uniformly  sampled,  discrete-time  sequences  have  been  reviewed.  These 
concepts  lead  to  the  generalization  of  the  pre-envelope,  involving  a  circular  frequency 
rotation  of  the  complex  envelope.  The  generalized  pre-envelope  along  with  the 
development  and  representation  of  the  modulo  2tc  process,  g[a,/J,  together  reveal 
themselves  to  be  fundamental  in  the  understanding  of  the  discrete-time  FM  (DTFM)  signal 
and  backward  difference  demodulation.  Subsequently,  a  method  of  enhancing  the 
performance  of  the  numerical  demodulation  process  using  frequency  feedback  is 
presented.  A  form  of  this  numerical  FM  with  feedback  (FMFB)  demodulator  is 
introduced  and  referred  to  as  the  reconstituted  numerical  FM  with  feedback  (RNFMFB) 
demodulator.  Simulation  results  are  presented  for  a  specific  implementation  of  this 
demodulator  operating  in  the  presence  of  noise.  Performance  enhancement  is 
demonstrated  for  large  Af/fm  ratios  at  low  sample  rates,  and  low  carrier-to-noise  ratios. 
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Regarding  the  backward  difference  numerical  FM  demodulation  process,  the 
presence  of  additive  white  Gaussian  noise  is  shown  to  create  an  additive  Middleton 
distributed  phase  noise  term,  along  with  randomly  occurring  +/-2tt  phase  cycle  changes  or 
phase  cycle-slips.  •  These  phase  cycle-slips  increase  in  rate  of  occurrence  with  an  increase 
in  the  modulation  index,  and  are  equivalent  to  noise-induced  instantaneous  frequency 
aliasing.  The  determination  of  the  rate  of  occurrence  of  these  slips  when  the  message 
waveform  is  known,  is  one  performance  prediction  measure  used  in  this  analysis. 
Experimental  results  obtained  via  simulation  of  the  noisy  DTFM  signal  agree  with  the 
analytically  predicted  phase  cycle-slip  occurrence  rates.  The  second  performance 
prediction  measure  is  the  estimate  of  phase  noise  variance.  Experimental  measurements  of 
this  variance  estimate  also  gave  results  in  agreement  with  the  true  variance  of  the 
Middleton  distribution.  These  performance  measures,  the  phase  cycle-slip  probability 
estimate  and  the  phase  noise  variance  estimate,  together  establish  a  baseline  performance 
summary  as  associated  with  backward  difference  numerical  FM  demodulation  techniques. 

The  sampling  rate  requirements  of  the  DTFM  signal  have  been  addressed.  Two 
parameters,  the  relative  message  sampling  rate,  y,  and  the  DTFM  modulation  index,  8,  are 
required  when  presenting  the  analysis  of  waveforms  composed  of  modulated  complex 
exponentials.  Proper  selection  of  these  parameters  ensures  proper  recovery  of  the  original 
message  from  the  complex  exponential  signal. 
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6.2  Areas  of  Future  Research 


In  this  research  we  have  intentionally  separated  the  phase  cycle-slip  error 
mechanism  from  the  Middleton  phase  noise  process,  to  elucidate  device  behavior  and  to 
lead  to  a  method  of  performance  enhancement.  With  careful  consideration,  these 
measures  can  be  combined  to  form  a  measure  of  mean  square  error,  leading  to  the  use  of 
optimal  receiver  design  and  analysis  techniques  [40-42],  Whether  working  in  the  phase  or 
instantaneous  frequency  domain,  it  may  further  be  appropriate  to  include 
pre-emphasis/de-emphasis  filtering  in  the  communication  system  model.  More  generally, 
wherever  additional  information  regarding  the  message  signal  is  available  a  priori,  it  is 
often  useful  to  modify  the  system  model  accordingly.  In  addition  to  presenting  results  in 
the  context  of  optimal  receiver  design  and  analysis  using  a  minimum  mean  square  error 
criterion,  the  non-ideal  effects  of  bandpass  filtering  on  complex  envelope  generation  can 
be  addressed.  Although  methods  of  generating  the  complex  generalized  pre-envelope 
from  the  received  band-limited  signal  have  been  presented,  of  interest  would  be  the  effects 
of  filtering  in  both  the  static  and  tracking  filter  cases.  These  filter  effects  are  particularly 
important  when  the  receiver  employs  frequency  feedback  techniques. 

Various  methods  exist  in  the  literature  which  help  to  alleviate  the  onset  of  phase 
cycle-slips,  which  tend  to  occur  (and  increase  rapidly)  as  the  input  carrier-to-noise  ratio 
decreases.  In  particular,  the  phase-lock  loop  device  has  been  a  popular  choice  for 
communications  engineers.  Also  of  interest,  however,  have  been  methods  of 
demodulation  which  employ  frequency  feedback  techniques.  Further  research  on  the 
development  and  application  of  the  RNFMFB  demodulator  can  be  performed. 
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Specifically,  the  inclusion  of  complex  bandpass  filtering  prior  to  discrimination,  and 
lowpass  filtering  in  the  feedback  path  are  expected  to  be  extremely  useful  in  various 
applications,  including  co-channel  interference  mitigation.  The  ability  to  separately 
control  the  distinct  function  of  these  filters  makes  the  FMFB  device  attractive  to 
communication  engineers.  In  addition,  methods  of  incorporating  additional  knowledge 
regarding  the  message  signal  into  the  phase  prediction  process  can  be  addressed. 

Other  areas  of  future  research  include  the  potential  application  of  these  results  to 
the  problem  of  phase  unwrapping  of  frequency  domain  responses  in  homomorphic 
processing  applications  [12,  43],  the  use  of  the  central  difference  derivative 
approximation,  data  signal  recovery  [44],  and  low  sample  rate  demodulation  of 
over-modulated  DTFM  signals. 
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APPENDIX  A 

The  Backward  Difference  Approximation  to  Time  Differentiation 
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A.l  Acquisition  of  a  Bandlimited  Signal  and  its  Corresponding  Derivative 

Consider  the  two  methods  of  acquiring  both  the  sampled  version  of  a  bandlimited 
signal,  x(t),  and  the  sampled  version  of  its  time  derivative,  x'(t),  as  shown  in  Figure  (A-l). 
The  sampling  time  interval,  Ts,  is  chosen  such  that  the  sampling  rate,  Fs  =  \ITS,  is  greater 
than  twice  the  highest  frequency  component  of  x(t).  As  indicated  in  the  figure  the  integer, 
n ,  is  an  index  which  identifies  a  particular  sampling  instant  via  the  product,  t  =  nTs. 

The  input,  x(t),  has  a  Fourier  transform  such  that 

\X(f)\  >  0,  0<\f\<W 

=  0,  else.  (A-l) 

Note  that  both  systems  are  considered  to  be  linear  and  time-invariant. 

A.2  The  Frequency  Response  of  the  Time  Derivative  Process 

For  System  A  the  ideal  time  derivative,  y(t),  is  known  to  have  the  Fourier 
transform  [29] 


Y(f)  =j2nf-X(f) ,  (A-2) 

which  is  also  zero  outside  the  interval  0  <  \f\  <  W,  since  X{f)  is  zero  outside  this  interval. 
This  system  requires  synchronized  sampling  of  both  the  input  signal,  and  its  corresponding 
derivative. 
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Figure  (A-l).  Two  methods  of  acquiring  a  bandlimited  signal,  x(t),  and  its  corresponding 
time  derivative,  a)  A  continuous-time  method  employing  an  analog  time  derivative  and 
synchronized  Analog  to  Discrete-time  Converter  (ADC)  pair,  b)  The  discrete-time 
backward  difference  approximation,  employing  a  single  ADC  and  difference  of 
consecutive  samples. 


We  can  represent  x(nTs)  as  the  superposition  of  infinitesimally  small  complex 
sinusoids  of  the  form  [12] 


±-X(eJQ)-eJQndn, 


as  in  the  Fourier  integral 


x(nTs)  =  ~  |  X{e^)e^nd£i  . 


(A-3) 


Note  that  Q  ranges  over  an  interval  of  length  271,  and  X{e?n)  determines  the  relative 
(complex)  amount  of  each  complex  sinusoidal  component.  Here,  Q  =  (oTs  radians,  and 
©  =  27 if  radians/second  1  .  Likewise,  y(nTs)  can  be  represented  in  the  above  manner. 

The  frequency  response  of  System  A  is  represented  as  HA(e’n),  where 

Y(eJQ)  =  HA(ejU)  •  X(e&)  .  (A-4) 

We  can  determine  HA(e?n)  by  considering  the  response  of  System  A  to  the  input 

x(t)  =  e*** ,  (A-5) 


where  ©  =  Q  /  Ts .  In  this  case  we  find  that 


=  {im1'11'  (  I MT: 

=  [/'S^  •  (A-6) 


Thus  HA(e,n)  can  be  identified  in  Equation  (A-6)  as 


1  To  be  consistent  with  the  usual  meaning  of  the  variable  ©,  the  variable  Q.  must 
have  units  of  radians.  (This  differs  from  the  usage  of  these  variables  in  reference  [12].) 
When  Ts  is  equal  to  1  second,  Q  and  ©  happen  to  have  the  same  value  and  it  is  tempting  to 
consider  them  as  having  the  same  units,  but  strictly  speaking,  the  units  remain  as  given. 
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HA(eJ(oTs)  =  HA(ejU) 


= jco  ,  -7i  < &TS  <  +7i  .  (A-7) 


Note  that  this  frequency  response  differs  from  that  implied  by  Equation  (A-2),  since 
HA(e?a)  is  implicitly  periodic  in  Q,  with  period  27t.  However,  since  x(t)  is  bandlimited  and 
properly  sampled,  the  effect  of  these  two  distinct  responses,  over  the  range  of  ©  of 
interest,  is  the  same. 

The  magnitude  frequency  response  of  System  A  is 


\HA(e^)\  =  jHA(e^)'HA(e-^) 

=  lool  ,  -n  <  (oTs  <  +7i  .  (A-8) 

The  phase  response  is  identified  as 


/lHA(e]aTs)  -  y  ,  -n  <(oTs  <0 
=  0  ,  &TS  =  0 
=  ^  ,  0  <  &TS  <  +7t 


(A-9) 
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A.3  The  Frequency  Response  of  the  Backward  Difference  Approximation 

In  a  similar  manner,  the  transfer  function  of  System  B  is  represented  as  HB{e?a), 

where 


Y(ejQ)  =  HB(ejQ )  •  X(Vn)  .  (A- 1 0) 

We  can  determine  HB(e?Q)  by  considering  the  response  of  System  B  to  the  input  identified 
by  Equation  (A-5).  For  System  B  we  find  that 


yinTs)  =  iWnr,)  -*([«  -  1  ]TS)} 
=  jr{ejQn  - 

=  jr{  \  -  e~jQ}  ■  eJQn  . 


(A-ll) 


From  Equation  (A-ll)  we  can  readily  identify  the  response  of  System  B  as 


=  y  ■  { 1  - e~J<oTs} ,  -7t  < (oTs  <  +7t  .  (A-12) 

*  s 

Equation  (A-12)  further  simplifies  to 

HB(et*T')  =  ys  sin(©n/2)  •  e*™7'*  ,  -it  < ®TS  <  +n  .  (A-13) 
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The  magnitude  frequency  response  of  system  B  is 


\HB(ejaT’)\  =  yJsm\®TJ2) 

=  jrj2~2cos(®Ts)  ,  -%<&TS  <  +71  .  (A-14) 


and  its  phase  response  is  identified  as 


ZHB(eJ<aTs)  =  (tc  —  (oTs)/2  -it,  -%< (oTs  <  0 
=  0  ,  (oTs  =  0 

=  (ti  -  <aTs)/ 2  ,  0  <  &TS  <  +7t  .  (A- 15) 

With  the  frequency  responses  of  the  ideal  time  differentiator  and  the  backward  difference 
having  been  determined,  it  is  possible  to  determine  how  closely  the  difference 
approximation  resembles  the  ideal  differentiator. 
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A.4  The  Equivalent  Distortion  Response 

We  can  consider  the  frequency  response  of  System  B  in  Figure  (A-l)  as  the 
product  of  the  ideal  response  of  System  A,  and  a  distortion  response,  HD{e?Q),  as  shown  in 
Figure  (A-2). 


HB{ejmTs)  =  HA{ejmlv)  •  HD(e^Ts)  . 


Figure  (A-2).  A  discrete-time  filter  which  has  the  same  frequency  response  as  System  B 
shown  in  Figure  (A-l). 


It  is  readily  seen  that  the  distortion  response  is 


HD(e»T') 


HB(e*°T') 
HA(e^)  ' 


(A- 16) 


From  Equation  (A-7)  with  j  =  e'7l/2 ,  along  with  Equation  (A- 13)  we  find  that 


HD(e*»T’) 


sin(cor,/2)  •  e1(n~aTs')'2 


=  sinc(ffl7V2)  •  e~^T^  ,  -k  <  & Ts  <  +k  ,  (A-l 7) 
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where  sinc(x)  is  defined  as  sin(x)  /  x  .  Since  sinc(x)  is  positive  over  the  interval 
-nil  <x<  +k/2  ,  the  magnitude  frequency  response  is  identified  as 

\Hj)(eJaT,)\  =  sinc(co7,.s/2) ,  -n  <  a>Ts  <  +n  ,  (A-18) 

and  is  periodic  in  Q  =  &TS ,  with  period  2n  .  The  phase  of  this  equivalent  distortion 
response  is 


ZHD{ejmTs)  =  ~g)Ts/2  ,  —71  <  oTs  <  +7r  .  (A-19) 

Equation  (A-19)  places  in  evidence  the  fact  that  the  backward  difference 
approximation  introduces  a  half-sample-interval  time  lag  relative  to  the  input  function  and 
relative  to  the  ideal  derivative  of  this  input.  Although  this  time  delay  is  not  a  distortion 
per  se,  knowledge  of  the  existence  of  this  delay  may  affect  design  considerations  in  some 
applications.  Thus  the  distortion  phase  response  is  linear,  and  not  generally  considered  to 
be  a  problem. 

The  magnitude  frequency  response  of  the  distortion  is,  in  fact,  not  particularly  bad 
either.  As  seen  in  Figure  (A-3),  the  magnitude  response  of  System  B  has  the  effect  of 
attenuating  higher  frequencies  as  compared  to  the  magnitude  response  of  the  ideal 
differentiator  of  System  A.  Figure  (A-4)  shows  the  magnitude  and  phase  response  of  the 
equivalent  distortion,  HD(e?°). 
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Figure  (A-3).  One  period  of  the  magnitude  frequency  responses,  Ts\HA(eJwTs)\  and 
rs\HB(e^)\  . 
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APPENDIX  B 


The  Backward  Difference  Phase  Estimate  Representation 


B.l  The  Effects  of  the  Presence  of  Phase  Noise 

We  are  interested  in  justifying  the  representation  of  our  noise-contaminated  phase 
estimate  as  the  sum  of  the  original  phase,  a (nTs),  the  Middleton  distributed  phase  noise, 
r\(nTs),  and  a  sheet  sequence  phase  error,  2%Re(nTs).  Also  of  interest  are  any  constraints 
or  properties  of  the  integer  sequence,  Re(nT),  that  may  arise. 

From  Figure  (4-2),  it  is  apparent  that  the  recovered  instantaneous  phase  estimate 

can  be  represented  as 


<j>,(w7V)  =  g[g[a(nTs)]  +  ri(«7;)] 

=  g[a(nTs )  +  r\(nTs)]  (B-l) 

(Here,  the  modulo  2n  arithmetic  property  presented  in  Equation  (3-49)  has  once  again 
been  utilized,  and  will  be  employed  without  reference  throughout  the  remainder  of  this 
appendix.)  By  definition  we  have  from  Equations  (4-20)  and  (4-21)  that 

$>r,)  =  Fs  ■  gElUnT.)  -?,([»-  1]  T,)]  (B-2) 


and 


UnTs)  =Ts-{£  %'AkTs)  +  M-n)} •  (B-3> 


Thus  from  Equations  (B-l)  and  (B-2)  we  can  write 
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Vd(nTs)  =  Fs  •  g[g[a(»7;>  +  y\(nTs)}  -  g[a([n  -  1]^)  +  r,([n  -  1]7,)]] 


=  Fs  ■  g[a'd(nTs)  +  t\ d(nTs)] ,  (B-4) 

where  a’d(nTs)  and  r\'d(nT)  are  as  defined  in  Equations  (4-38)  and  (4-39)  respectively.  In 
the  absence  of  noise,  to  prevent  instantaneous  frequency  aliasing,  we  have  from  Chapter  3 
the  requirement 

-7t  <  a'd(nTs )  <  +7t  .  (B-5) 

The  phase  noise,  ii,  has  been  defined  such  that 

— 7t  <  r\(nTs)  <  +7t , 

which  results  in  the  property 

-271  <  T\d(nTs)  <  + 2 7t  .  (B-6) 

Combining  Equations  (B-5)  and  (B-6)  we  have  that 

-37t  <  ad(nTs )  +  r\'d(nTs)  <  +3%  .  (B-7) 

From  Equation  (B-7)  we  can  write 

ad(nTs)  +  r\d(nTs)  =  g[a/c/(«7’s)  +  r,'  (nTs)]  -  2%R'e(nTs)  ,  -1  <  R'e(nTs)  <  +1 .  (B-8) 
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Note  that  for  any  particular  time,  t  =  nTs,  we  have  that  R'e(nTs )  e  {-1,0, 1 } . 

If  we  had  available  to  us  the  quantity  a'd(nT)  +  r\'d(nT),  along  with  initial  phase 
conditions,  we  could  form  the  accumulation 


$goai(nTs)  =  E  {a ’d(kTs)  + 1) 'd{kTs)}  +  {a (~TS)  +  r\(-Ts)} 

k= 0 

=  a(nTs)  +  i\(nTs) .  (B-9) 


In  lieu  of  Equation  (B-9),  we  use  what  is  available  and  form 

fa(nTs)  =  E  g[a'd(kTs)  +  n 'd(kTs)\  +  $,(-/’,) ,  (B-lOa) 


where 


Combining  Equations  (B-8)  through  (B-10)  we  have  that 


(B-lOb) 


$goai(nTs)  =  $d(nTs)-2n  i  R'eQcTs)  ,  -1  <R'e{kTs)<+ 1  .  (B-ll) 


Defining  Re(nT)  such  that  R'e{nTs )  =  Re(nT ')  - ./?,([«-!] 7^),  we  arrive  at 


$goai(nTs)  =  -  2KRe(nTs)  ,  (B-12) 


where  the  integer  sequence  Re{nTs )  has  the  property  that 
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-1  <Re(nTs)-Re(\n-X\Ts)<+\  . 


(B-13) 


From  Equations  (B-9)  and  (B-12)  we  can  write 

§d(nTs)  =  a (nTs)  +  r|(«7;)  +  2%Re(nTs )  .  (B-14) 

It  is  also  readily  apparent  from  Equations  (B-l)  and  (B-14)  that 

UnTs)  =  g$d(tiTs)\  .  (B-l 5) 

This  fact  is  utilized  in  Section  4.1.3  to  help  establish  a  pair  of  experimental  statistics  as 
measures  of  backward  difference  phase  estimate  performance  prediction. 

Thus  our  noise-contaminated  backward  difference  phase  estimate  is  represented  as 
the  sum  of  the  original  phase,  a (nT),  the  Middleton  distributed  phase  noise,  r|(«r),  and  a 
sheet  sequence  phase  error,  2nRe(nTs).  The  sheet  sequence  error,  Re(nTs),  has  the 
property  that  it  either  remains  unchanged  between  consecutive  time  samples,  or  changes 
by  +/- 1 .  Therefore  we  refer  to  the  sequence 

2nRcs(nTs)  =  2k. R'e(nTs)  =  2k  •  {Re(nTs)  - Re([n -  l]r,)}  (B-16) 

as  a  phase  cycle-s/zp  sequence,  since  it  indicates  the  occurrences  of  +/-2k  phase  changes  at 
any  particular  time,  t  =  nTs. 
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